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OVERVIEW  OF  RESEARCH 


1.  Introduction 

The  research  conducted  under  this  contact  has  been  concerned  principally  with  the 
development  of  numerical  models  of  the  Earth's  low-latitude  boundary  layer  (LLBL),  a  thin  layer 
of  antisunward  flowing  plasma,  located  immediately  Earthward  of  the  equatorial  magnetopause 
current  layer  which  marks  the  outer  boundary  of  the  Earth's  magnetic  Held.  The  work  has  been 
concentrated  in  five  areas: 

1.  Development  of  a  self-consistent  steady-state  numerical  model  of  the  equatorial 
portion  of  the  LLBL  on  closed  field  lines,  including  coupling  to  the  ionosphere  via 
field-aligned  currents. 

2.  Development  of  a  self-consistent  numerical  model  of  the  force-free  boundary 
layer  that  provides  the  link  between  the  equatorial  LLBL  and  the  dayside  auroral 
ionosphere. 

3.  Examination,  by  use  of  numerical  simulation,  of  the  stability  of  laminar  flow  in 
the  equatorial  LLBL  in  the  presence  of  coupling  to  the  ionosphere  and  associated 
nonuniform  bending  of  the  magnetic  field  lines  in  the  LLBL. 

4.  Examination  of  resistive  tearing-mode  instability  in  a  current  sheet  with 
equilibrium  viscous  stagnation-point  flow. 

5.  Examination  of  magnetic-field  maxima  observed  in  the  low-latitude  boundary 
layer. 

In  the  following  sections  of  the  report,  a  brief  summary  and  discussion  of  the  results 
obtained  in  each  area  is  provided.  Details  of  the  research  are  provided  in  five  papers,  three  that 
have  been  published  and  two  more  that  have  been  submitted  for  publication.  Also  provided  are 
appropriate  extracts  from  E.  Drakou's  Ph.D.  thesis  which  is  concerned  with  the  LLBL  model  and 
was  developed  with  principal  support  from  the  present  contract.  All  of  these  documents  are 
appended  to,  and  form  an  integral  part  of,  this  final  report.  The  research  described  in  the  five 
papers  has  also  received  partial  support  from  other  sources,  as  indicated  in  their  acknowledgment 
sections. 
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2.  Numerical  Model  of  the  Equatorial  Low-Latitude  Boundary  Layer  on  Closed 
Field  Lines 

The  low-latitude  boundary  layer  (LLBL)  is  a  narrow  region,  located  in  the  low  latitude 
region  immediately  inside  the  outer  boundary  of  the  magnetosphere,  die  magnetopause.  The  LLBL 
contains  plasma,  mostly  of  magnetosheath  origin,  that  flows  along  the  layer  in  the  antisolar 
direction  at  a  speed  comparable  to  the  magnetosheath  flow  speed.  This  plasma  flow  is  at  an  angle 
—  in  the  simplest  model  at  a  90*  angle  —  to  the  geomagnetic  field  in  the  vicinity  of  the  equatorial 
plane  and  thus  it  has  an  associated  convection  electric  field,  Ee,  which  is  projected,  in  pan  at  least, 
into  the  ionosphere  at  the  footpoints  of  the  geomagnetic  field  lines  threading  the  LLBL.  This 
impressed  electric  field,  Ei,  drives  a  horizontal  Pedersen  current,  Ji,  in  the  ionosphere;  the 
divergence  in  E;  gives  rise  to  a  corresponding  divergence  in  this  horizontal  current,  i.e.,  it  gives 
rise  to  a  corresponding  magnetic-field-aligned  current  into  or  out  of  the  ionosphere.  This  field- 
aligned  current  connects  the  ionosphere  to  the  LLBL,  where  it  is  again  deflected  to  form  a  current 
Je  that  flows  across  the  equatorial  geomagnetic  field.  In  the  ionosphere,  the  product  Ej-Jj  >  0, 
whereas  in  the  LLBL  the  product  Ee-  Jc  <  0;  thus  the  former  region  acts  as  an  electrical  load  and 
the  latter  region  as  an  electrical  generator,  connected  to  the  ionospheric  load  via  the  field-aligned 
currents.  In  the  simplest  conceptual  model,  the  projection  of  the  equatorial  electric  field  into  the 
ionosphere  occurs  by  assuming  the  geomagnetic  field  lines  to  be  equipotentials.  In  more  realistic 
modeling,  a  potential  drop,  AO|j,  along  the  field  lines  is  incorporated  via  a  field-aligned 
conductance  K,  so  that  Jg  =  KAOjj.  In  the  post-noon  LLBL,  the  field-aligned  current,  J||,  flows 
out  of  the  post-noon  ionosphere  so  that  the  potential  drop  A<X>g  can  accelerate  electrons  precipitating 
into  the  ionosphere  to  energies  comparable  to  those  needed  to  explain  auroral  emissions.  On  the 
pre-noon  side,  a  potential  drop  AO||  will  accelerate  electrons  upwards  and  ions  downwards  instead. 
A  schematic  drawing  of  the  dawnside  LLBL  configuration  and  its  coupling  to  the  ionosphere  is 
shown  in  Figure  1.  The  equatorial  portion  of  the  LLBL,  in  which  plasma  inertia,  pressure  and 
viscosity  are  important,  is  located  in  the  region  Izl  $  H.  The  force-free  coupling  region,  where  the 
currents  are  entirely  field  aligned,  connects  to  the  LLBL  at  Izl  =  H  and  then  extends  along  magnetic 
field  lines  into  the  northern  and  southern  ionospheres.  Note  that  the  main  field-aligned  current 
associated  with  the  LLBL  provides  the  portion  of  the  so-called  Region  1  current  that  is  observed  to 
flow  into  the  pre-noon  (8  -  12  LT,  say)  and  out  of  the  post-noon  (12  -  16  LT,  say)  sides  of  the 
dayside  auroral  oval.  Any  field-aligned  return  current,  at  the  outer  edge  of  the  LLBL,  i.e.,  at  the 
magnetopause  itself,  could  correspond  to  the  so-called  NBZ  currents,  observed  at  low  altitudes 
during  conditions  of  northward  interplanetary  magnetic  field  (IMF).  Note  also  that  the  actual  local 
time  extent  of  the  LLBL  projection  into  the  ionosphere  is  unknown  a  priori  and  must  be  calculated 
in  a  self-consistent  manner  from  the  currents  in  the  LLBL  and  in  the  coupling  region. 
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Figure  1.  View  from  the  sun  of  the  dawnside  low-latitude  boundary  layer  and  its  coupling 
to  the  dayside  auroral  oval.  Coordinates  (x,y,z)  are  the  usual  GSM  coordinates. 

Finally,  it  is  emphasized  that  the  boundary  layer  is  assumed  to  be  located  on  closed 
geomagnetic  field  lines,  i.c.,  field  lines  that  have  both  ends  rooted  in  the  Earth.  There  is 
observational  evidence  that,  at  least  during  periods  of  southward  interplanetary  magnetic  field 
(IMF),  the  portion  of  the  LLBL  immediately  adjoining  the  subsolar  magnetopause  is  on  open  field 
lines  as  a  result  of  reconnection.  The  model  developed  under  this  contract  does  not  apply  to  this 
portion  but  it  would  apply  to  any  remaining  part  of  the  subsolar  LLBL,  which  is  on  closed  field 
lines:  in  such  a  situation,  the  outer  edge  of  the  LLBL  described  by  the  model  would  be  located  at 
the  first  closed  field  line  (i.e.,  the  inner  separatrix  of  the  reconnection  configuration)  rather  than  at 
the  magnetopause  itself.  For  northward  IMF  it  is  expected  that  most  or  ail  of  the  LLBL  will  be  on 
closed  field  lines;  this  may  also  be  the  case  on  the  magnetospheric  flanks,  regardless  of  IMF 
direction. 

The  model  of  the  equatorial  portion  of  the  LLBL  that  has  been  developed  under  the  present 
contract  is  an  outgrowth  of  earlier  one-dimensional  analytic  descriptions  of  the  LLBL  on  closed 
field  lines,  developed  by  Sonnerup  [1980],  Lotko  et  al.  [1987]  and,  in  particular,  Phan  et  al. 
[1989].  The  new  model  is  numerical  rather  than  analytic.  It  represents  the  simplest  possible 
extension  of  the  Phan  et  al.  model  to  include  evolution  of  boundary  layer  structure  and  thickness  in 
the  main  flow  direction  (-x)  along  the  magnetopause.  From  a  computational  viewpoint  the  model 
is  two  dimensional  in  the  sense  that  the  equations  describe  behavior  in  the  equatorial  (xy)  plane. 
However,  variations  with  the  third  coordinate  (z)  are  described  to  lowest  order,  including  a  self- 
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consistently  calculated  bending  of  the  magnetic  Held  lines  in  the  xz  plane  as  well  as  diamagnetic 
changes  of  the  field  component  Bz.  This  field  deformation  results  from  currents  flowing  in  the  y 
direction  across  the  boundary  layer  and  being  gradually  deflected  to  field-aligned  currents  that  flow 
into  the  ionosphere;  diamagnetic  currents  also  flow  in  the  x  direction.  The  height,  H,  (in  the  z 
direction)  of  the  boundary  layer  above  and  below  the  equatorial  plane  should  also  be  allowed  to  be 
a  function  of  x  and  y,  rather  than  being  constant  as  in  the  Phan  et  al.  model,  although  this  feature 
has  not  yet  been  incorporated.  The  condition  that  should  be  used  to  determine  H  is  that,  at  z  =  ±H, 
the  boundary-layer  plasma  pressure  has  dropped  off  to  its  ambient  magnetospheric  value. 

The  model,  which  is  described  in  detail  in  Appendix  1,  is  based  on  four  important 
simplifications  of  the  full  3D  MHD  equations;  (1)  steady  flow;  (2)  the  boundary  layer 
approximation  d/dy»d/dx  which  leads  to  kinematic  treatment  of  the  transverse  (y)  motion  and  to 
neglect  of  magnetopause  curvature  effects;  (3)  kinematic  treatment  of  the  plasma  motion  in  the  ±z 
direction  towards  or  away  from  the  equatorial  plane;  (4)  division  of  the  model  into  three  modules, 
namely  the  equatorial  LLBL  module  (the  current  generator),  the  force-free  coupling  module,  and 
the  ionospheric  module  (the  resistive  load).  At  present,  each  module  contains  only  the  simplest 
description  of  the  most  important  physics.  As  mentioned  already,  self-consistent  variations  of 
pressure,  density  and  magnetic  field  with  die  coordinate  z  are  included  but  the  main  computation  is 
in  effect  two-dimensional,  dealing  only  with  quantities  evaluated  at  z=0.  The  equations  governing 


the  complete  model  are 

pOVO'Vvx  =  -  dPoo(x)/dx  +  BzoBxi/poHr  +  (9/3y)(r|8vx/3y),  (1) 

PO  +  Bi/2W)  =  Poo(x),  (2) 

H(x,y)/Hr  =  { 2po[po(x,y)  -  PooWJ/B^fx.y) } 1/2,  (3) 

V(poHvo)  =  0,  V(Bzovo)  =  0,  vo-V(po/po"0  =  0.  (4) 


The  quantities  po,  po  and  Bzo,  along  with  the  velocity  vo  =  (vx,vy,0)  are  evaluated  in  the  equatorial 
plane;  p»(x)  is  the  magnetospheric  plasma  pressure,  T|  is  the  viscosity  (which  may  be  of  either 
microscopic  or  turbulent  origin),  and  Bx(x,y,z)  =  Bxi(x,y)(z/Hr),  Hr  being  a  constant  reference 
value  of  H.  This  z  dependence  of  Bx  leads  to  approximately  parabolic  field  lines  in  the  model, 
with  field-line  curvature  that  varies  with  the  coordinate  y,  being  a  maximum  at  the  magnetopause 
and  then  decreasing  as  one  moves  further  Earthward.  In  order,  the  above  equations  express: 
momentum/force  balance  in  the  x,  y,  and  z  directions  (the  expression  for  H(x,y)  given  in  the  third 
equation  is  derived  from  p  +  Bx2/2po  =  P0(x,y));  mass  conservation;  flux  conservation;  and 
isentropic  compression/expansion,  respectively.  To  these  equations  are  added  jxB  =  0  in  the 
coupling  module,  and  the  ionospheric  laws  jm  =  K(d>e  -  d>j)  s  KAOh,  and  d/9y(XpdC>j/9y)  =  -  jm. 
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where  <X>e(x,y,H)  and  d>i(xj,yi)  are  the  potential  distributions  in  the  LLBL  and  in  the  ionosphere, 
respectively.  In  the  simplified  version  of  the  model  produced  under  this  contract,  the  LLBL  height 
H  is  constant,  the  conductance,  K,  along  field  lines  is  infinite,  and  the  ionospheric  conductivity, 
£p,  is  constant.  Relations  between  the  coordinates  (x,y)  and  (xi,yi),  expressed  via  mapping 
factors,  can  be  calculated  self -consistently  in  the  coupling  region,  although  in  die  current  version  of 
the  model  the  mapping  factors  are  taken  to  be  constant.  In  other  words,  the  magnetic  field 
configuration  in  the  LLBL  itself,  as  well  as  in  the  force-free  coupling  region,  would  ultimately  be 
computed  self  consistently  (albeit  in  the  boundary  layer  approximation).  Thus,  a  complete  version 
of  the  model  would  provide  an  accurate  mapping  along  B  from  the  ionosphere  to  the  equatorial 
plane.  The  field  from  a  realistic  model  of  the  inner  magnetospheric  B  field  (e.g.,  one  of  the 
Tsyganenko  models)  would  then  be  used  as  a  boundary  condition  at  the  magnetospheric  edge  of 
the  LLBL  computational  box.  The  mapping  in  the  coupling  region  and  the  interface  with  a  suitable 
magnetospheric  field  model  are  described  in  Section  3  and  in  Appendix  S. 

The  boundary  layer  equations  are  parabolic  and  are  therefore  integrated  by  use  of  a 
marching  procedure  (in  the  -x  direction)  that  allows  one  to  follow  the  development  of  the  layer  for 
as  large  distances  in  the  flow  direction  as  available  computer  resources  permit.  Since  the 
downstream  state  obtained  from  one  run  can  be  used  as  the  upstream  state  for  a  second  run,  etc., 
the  boundary  layer  evolution  can  in  principle  be  followed  to  arbitrarily  large  distances  along  the 
magnetospheric  flank.  One  of  the  difficulties  associated  with  the  numerical  marching  procedure  is 
that  it  cannot  handle  reversals  in  the  main  flow  direction,  i.e.,  in  the  velocity  component  vx.  This 
difficulty  is  associated  with  the  inertia  term  pvxdvx/dx,  which  does  not  reverse  sign  when  vx 
reverses  sign.  The  situation  is  similar  to  the  integration  of  a  diffusion  equation  backwards  in  time 
which  is  numerically  unstable.  This  problem  has  been  overcome  in  die  present  code  by  neglecting 
inertia  near  the  flow  reversal  and  in  the  entire  region  of  slow  sunward  flow  on  the  magnetospheric 
side  of  the  LLBL.  Furthermore,  in  the  present  version  of  the  code,  the  velocity  profile,  vx(x,y),  in 
the  slow-flow  region  is  obtained  analytically  by  assuming  that  the  plasma  properties  and  the  field 
component  Bz  have  reached  their  asymptotic  magnetospheric  values  in  this  region.  In  other  words, 
they  depend  on  x  but  not  on  y.  This  procedure  decreases  the  size,  along  y,  of  the  computational 
domain,  thus  providing  for  substantial  computational  economy.  Benchmark  tests  (Appendix  4)  in 
which  solutions  from  the  program  are  checked  against  certain  exact  self-similar  solutions  indicate 
that  the  procedure  adopted  provides  satisfactory  accuracy. 

Details  of  the  numerical  procedure,  along  with  benchmark  tests  and  a  discussion  of  various 
generalizations  of  the  model,  are  given  in  Appendices  3, 4  and  2,  respectively. 

Sample  results  from  a  test  run  of  the  numerical  code  (with  H  =  const,  Ad>n  =  0  and  constant 
mapping  factors,  which  represents  its  current  status)  are  shown  in  Figure  2.  A  complete 
presentation  of  results  to  date  is  given  in  the  article  by  Drakou  et  al.  [1994]  which  is  reproduced  in 
Appendix  1. 
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Figure  2.  Results  of  trial  run,  using  the  present  boundary  layer  model. 

The  boundary  conditions  used  in  Figure  2  include  an  accelerating  tailward  plasma  flow  speed, 
vx(x,0),  at  the  magnetopause,  a  small  constant  sunward  convection  speed,  vxoo,  in  the 
magnetosphere,  where  plasma  number  density,  n«,  pressure,  p«,  and  temperature,  Too,  as  well  as 
Bz  Field  component,  BZOo,  all  decrease  with  increasing  distance,  -x,  down  the  tail.  Viscosity, 
assumed  in  this  run  to  be  proportional  to  p/B,  is  included.  Entrainment  by  the  LLBL  of 
magnetospheric  plasma  is  illustrated  in  the  top  middle  panel  where  the  vy  velocity  component  has 
been  artificially .  ugmented  in  order  to  make  the  effect  visible.  Viscous  widening  of  the  velocity 
profile  vx(x,y)  is  seen  in  the  middle  left  panel,  whereas,  n,  T,  and  Bz  (bottom  panels)  have  profile 
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widths  that  decrease  as  one  moves  downstream.  The  field-aligned  current  shows  considerable 
structure  (middle  right  panel)  and  substantial  evolution  (top  right  panel)  in  the  flow  direction.  Note 
in  particular  that  a  maximum  in  jz  occurs  at  a  particular  x  value,  i.e.,  at  a  particular  local  time;  such 
a  maximum  has  indeed  been  reported  in  the  observational  study  by  Iijima  and  Potemra  [1978]. 
These  results  are  shown  here  for  illustrative  purposes  only;  their  details  are  expected  to  change 
considerably  when  additional  features  such  as  field-aligned  potential  drops  are  incorporated. 

3.  Numerical  Model  of  the  Force- Free  Coupling  Module 

In  the  narrow  region  that  connects  the  equatorial  LLBL  to  the  ionosphere,  the  plasma 
density,  pressure,  and  velocity  are  sufficiently  small  and  the  magnetic  field  intensity  sufficiently 
high  so  that  the  electric  currents  connecting  the  LLBL  to  the  ionosphere  must  be  Held-aligned  to  a 
high  degree  of  accuracy.  Thus  the  basic  governing  equations  in  this  region  are 

jxB  =0,  j  =  (l/po)VxB,  V  B  =  0.  (5) 

These  equations  lead  to 

Poj  =  a(x,y,z)B  =  VxB  (6) 

where  a  is  a  proportionality  factor  which  has  a  constant  value  along  any  magnetic  field  line  so  that 

B  Va  *  0  (7) 

The  latter  two  relations  have  been  simplified  by  use  of  the  boundary  layer  approximation 
(described  in  Section  2)  and  a  computer  code  has  been  developed  for  generating  self-consistent 
current  and  magnetic  field  configurations  in  the  connection  layer.  Boundary  conditions  for  the 
code  consist  of  specification  of  the  magnetic  field  at  the  Earthward  edge  of  the  layer  by  use  of  an 
empirical  magnetospheric  field  model  such  as  one  of  the  Tsyganenko  models  and  specification  of  a 
field-aligned  current  distribution  and  magnetic  field  at  the  interface  to  the  LLBL  at  Izl  =  H 
(alternatively,  these  quantities  can  be  specified  in  the  ionosphere).  Ultimately,  this  code  should  be 
appended  to  the  main  LLBL  code  described  in  Section  2  but  this  step  has  not  yet  been  taken. 
When  the  two  codes  have  been  combined,  the  resulting  code  will  be  able  to  give  accurate 
information  about  the  magnetic  field  mapping  in  the  regions  immediately  Earthward  of  the 
magnetopause  where  existing  field  models  such  as  the  various  Tsyganenko  models  fail  to  provide 
reliable  information.  However,  the  utility  of  the  force-free  boundary  layer  module  is  not  restricted 
to  the  LLBL:  it  has  independent  applications  to  any  field-aligned  current  sheets  in  the 
magnetosphere. 

A  detailed  description  of  the  force-free  boundary  layer  model  is  provided  in  Appendix  5, 
along  with  a  benchmark  test  and  a  simple  application  using  measured  Region  1  currents  above  the 
ionosphere  [Iijima  and  Potemra,  1978]  in  conjunction  with  a  dipolar  magnetospheric  field.  These 
results  suggest  typical  azimuthal  deflections  of  dayside  magnetic  field  lines  in  the  range  21-26°. 
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4.  Kelvin-Helmholtz  Instability  in  the  Low-Latitude  Boundary  Layer 

The  stability  of  the  plasma  flow  in  the  LLBL  is  of  importance  for  the  construction  of  any 
realistic  model  of  this  layer.  Stable  behavior  implies  that  any  transport  of  momentum  or  mass 
across  the  layer  is  caused  by  microscopic  plasma  processes  while  unstable  behavior  implies 
transport  by  eddy  viscosity  and  eddy  diffusivity.  To  date,  studies  of  the  Kelvin-Helmholtz  (KH) 
instability  relevant  to  the  LLBL  have  assumed  the  unperturbed  magnetic  field  lines  to  be  straight. 
However,  in  a  realistic  model  of  this  layer,  the  currents  flowing  across  it  and  being  gradually 
deflected  into  field-aligned  currents,  as  the  Earthward  edge  of  the  layer  is  approached,  produce 
bending  of  the  Held  lines  into  approximately  parabolic  shapes  with  vertices  pointing  in  the 
(antisolar)  flow  direction  (see  Figure  5  of  Appendix  1).  Parabolas  close  to  the  magnetopause  have 
larger  curvature  than  those  close  to  the  magnetospheric  edge  of  the  LLBL  so  that  one  may  expect 
interchange  motions  to  be  impeded. 

We  have  carried  out  a  fully  three-dimensional  simulation  of  a  velocity  shear  layer  of  the 
hyperbolic  tangent  type  in  the  presence  of  coupling  to  the  northern  and  southern  ionospheres, 
represented  in  the  simulation  model  by  two  parallel  electrically  conducting  plates.  The  current 
system  produced  by  this  coupling  leads  to  parabolic  field  lines  of  opposite  curvature  on  the  two 
sides  of  the  shear  layer. 

The  results  of  the  simulation  are  reported  in  detail  in  Appendix  6.  In  brief,  it  is  found  that 
the  field  line  curvature,  if  substantial,  may  severely  suppress  the  KH  instability.  When  the 
curvature  is  less  strong,  as  in  the  dayside  LLBL,  the  instability  proceeds  albeit  at  a  somewhat 
reduced  rate  compared  to  the  case  of  straight  field  lines  transverse  to  the  flow,  and  leads  to  the 
formation  of  three-dimensional  vortex/current  structures  that  may  be  related  to  observed  auroral 
bright  spots. 

5.  Resistive  Tearing  Mode 

The  behavior  of  the  resistive  tearing  mode  near  the  subsolar  magnetopause  has  been 
investigated  with  the  objective  of  finding  out  how  the  presence  of  stagnation  point  flow  and 
viscosity  influence  this  instability.  The  unperturbed  equilibrium  is  an  exact  solution  of  the 
incompressible  MHD  equations,  including  resistivity  and  viscosity  so  that  the  stability  properties 
can  be  investigated  even  for  small  magnetic  and  viscous  Reynolds  numbers  where  the  traditional 
tearing  mode  analysis  is  invalid.  The  results  indicate  stability  of  the  tearing  mode  for  magnetic 
Lundquist  numbers,  S,  (based  on  the  Alfv6n  speed)  less  than  12.25,  regardless  of  viscous 
Reynolds  number,  Re,  and  for  S  £  18.3  when  Re  =  0.  For  large  values  of  S  and  Re,  the  classical 
asymptotic  growth  rates  of  Furth,  Killeen  and  Rosenbluth  are  recovered.  It  is  also  found  that  the 
stagnation  point  flow  stabilizes  long  wavelength  perturbations  so  that  the  tearing  mode  has  a  cut¬ 
off  at  small  as  well  as  at  large  wave  numbers.  The  main  effect  of  viscosity  is  to  reduce  the  growth 
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rate  of  the  instability,  in  particular  at  short  wavelengths.  The  stabilization  at  low  magnetic 
Reynolds  numbers  is  relevant  to  the  subsolar  magnetopause  where  estimates  of  S  are  in  the  range  2 
<  S  <  100  [Lee  and  Fu,  1986]. 

6.  Magnetic  Field  Maxima  in  the  LLBL 

It  has  been  known  for  a  long  time  that  the  magnetic  field  often  exhibits  a  maximum 
immediately  Earthward  of  the  magnetopause  instead  of  the  minimum  one  would  expect  from  the 
diamagnetic  effect  of  the  dense  plasma  in  the  LLBL.  We  have  investigated  this  effect  by  use  of 
plasma  and  magnetic  field  data  from  the  AMPTE/IRM  spacecraft  and  find  two  fundamentally 
different  causes  for  the  excess  field:  (i)  a  depression  within  the  LLBL  of  the  density  of  medium 
energy  ions  of  magnetospheric  origin  and  (ii)  field  curvature  effects  associated  with  undulations  of 
the  magnetopause  itself.  When  case  (i)  is  at  hand,  medium-energy  electron  fluxes  are  also 
depressed,  suggesting  that  the  field  lines  in  the  LLBL  may  have  been  open.  However,  this  is  not 
the  only  possible  explanation  for  the  absence  of  energetic  particles. 
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Remits  ate  presented  from  a  steady  state  numerical  model  of  the  low-latitude  boundary  layer  (LLBL)  on  closed 
field  lines  and  its  coupling  to  the  dayside  auroral  ionosphere.  In  the  model  the  boundary  layer  approximation  is 
used,  the  result  being  that  inertial  forces  are  taken  into  account  only  in  the  main  flow  direction  (— r)  where  they 
are  balanced  by  pressure  fortes,  j  x  B  forces,  and  viscous  forces.  Motion  in  the  transverse  directions  (y  and  a)  is 
beared  tcingmaHcally,  the  force  balances  in  them  two  directions  being  purely  static.  Computatiooally.  the  model 
is  two  dimensional,  describing  the  motion  of  plasma  and  frozen-in  magnetic  field  in  the  equatorial  (xy)  plane 
but  allowing  for  lowest-order  polynomial  variation  of  some  quantities  with  the  coordinate  (a)  perpendicular  to 
that  plane.  The  plasma  expands  and  compresses  isentropically;  the  magnetic  field  is  calculated  self-coosisrently. 
which  leads  to  approximately  parabolic  field  line  shape  in  planes  parallel  to  the  magnetopause  (the  xx  plane), 
with  maximum  field  curvature  near  the  magnetopause  edge  of  the  LLBL  Coupling  to  the  ionosphere  via  region  1 
field-aligned  currents  is  included.  The  effects  of  the  ionosphere  are  represented  by  two  parallel  resistive  plates  at 
fixed  height  above  and  below  the  equatorial  plane.  The  model  can  be  used  to  investigate  the  influence  of  various 
physical  parameters,  for  example,  viscous  and  magnetic  Reynolds  numbers,  and  of  boundary  conditions  at  the 
magnetopause  and  in  the  magnetosphere  on  the  LLBL  development  in  the  -x  direction.  Special  attention  is 
given  to  viscous  effects  which,  under  suitable  circumstances,  lead  to  a  region  I  current  that  first  increases  and 
then  reduces  with  increasing  longitude  away  from  local  noon.  Asymptotic  matching  of  the  antisunwatd  motion 
of  the  cool  LLBL  plasma  to  sunward  convection  of  hot  plasma  in  the  magnetosphere  is  illustrated  along  with  the 
entrainment  of  magnetospheric  plasma  by  the  antisunward  LLBL  flow. 


Introduction 

The  low-latitude  boundary  layer  (LLBL)  is  a  narrow  region  of 
tailwatd  flowing  plasma  located  in  the  magnetosphere,  immedi¬ 
ately  inside  the  magnetopause  current  sheet,  at  low  geomagnetic 
latitudes.  It  was  first  observed  by  Hones  et  aL  [1972]  and  Miasofu 
etal  [1973]  along  the  flanks  ofthe  geomagnetic  tail  Since  then,  a 
mtmberof  authors  have  discussed  observational  data  and  theoretical 
models  of  the  LLBL  From  these  studies  we  know  the  following: 
the  LLBL  is  intermittently  present  at  almost  all  local  times  along 
the  entire  dayside  portion  ofthe  magnetopause.  It  was  also  recently 
found  to  be  extended  to  evening  local  times  [Woch  etal,  1993].  Its 
thickness  appears  to  increase,  on  average,  with  increasing  distance 
from  the  subsolar  point  [Eastman  and  Hoots,  1979],  although  the 
growth  may  be  slow  beyond  the  dawn-dusk  meridian  plane.  A  typ¬ 
ical  boundary  tayer  thickness  at  that  location  is  0.5-1  Re  [Eastman 
and  Hones,  1979;  Sckopke  et  aL,  1981].  When  the  interplane¬ 
tary  magnetic  field  is  north  ward,  the  LLBL  appears  thicker  and  is 
thought  to  be  located  on  closed  field  lines,  that  is,  on  magnetic  field 
lines  that  have  both  feet  in  the  ionosphere  [HaerendeletaL,  1978; 
Williams  etaL,  19$5;  Mitchell  et  aL,  1987;  Trover  etaL,  1991].  It 
is  this  kind  of  LLBL  we  are  modeling  in  the  present  paper.  During 
periods  of  southward  interplanetary  magnetic  field  the  situation  is 
less  clear  reconnection  may  dominate  in  which  case  the  LLBL  on 
closed  field  lines  may  be  mostly  absent  The  experimental  evidence 
in  this  case  is  not  conclusive,  bu|  it  indicates  that  the  layer  is  par- 


1  Now  at  Hetzbeig  Institute  of  Astrophysics.  National  Research  Council 
of  Canada.  Ottawa.  Ontario. 

Copyright  1994  by  As  American  Geophysical  Union. 


Ptoer  number  93JA02094. 
0l4t4227/94/93JA42094$3.00 


TtwUjS^Oovwnment  Is  authorised  to  reproduca  and  sail  this  repot 
"""d"*'®"  by  others  must  be  obtained  fro 

wuynpu  owner. 


13 


daily  on  open  and  partially  on  closed  field  lines:  our  model  applies 
only  to  the  part  on  closed  field  lines.  As  one  moves  inward  from 
the  magnetopause  across  the  LLBL  the  density  falls  from  a  high 
magnetosheath  like  value  to  a  low  magnetospheric  value,  while  the 
magnetic  field  usually  rises  slightly  (although  a  field  enhancement 
is  sometimes  teen  instead  [Sonnerup  et  aL.  1992])  and  the  plasma 
pressure  falls  accordingly  [Sckopke  et  aL ,  1981].  The  dense  plasma 
in  the  LLBL  appears  to  be  mainly  of  magnetotheath  origin.  How 
this  plasma  enters  the  LLBL  is  not  yet  clear,  h  may  leak  diffusively 
over  portions  of  the  magnetopause  surface  [Tsurutani  and  Thome, 
1982]  or  enter  onto  closed  field  lines  in  the  cusps  [Pluchmann  et 
aL,  1976]  or  at  the  edges  of  a  dayside  reconnection  segment  of 
the  magnetopause.  It  is  also  possible  that,  during  northward  IMF, 
magnetosheath  flux  tubes  reconnect  in  the  norib  and  south  beyond 
the  cusps  and  then  shorten  and  reorient  themselves  to  eventually 
become  assimilated  with  other  magnetospheric  field  lines  [Dmtgey, 
1963;  Cowley  and  Owen,  1989;  Song  and  RnsseU,  1992].  As  noted 
by  Mitchell  etaL  [\9V1],SckopkeetaL  [1981J,  Eastman  and  Hones 
[1979],  and  others,  the  cool,  dense  magnetosbeath  like  plasma  in 
the  LLBL  is  usually  mixed  to  some  extern  with  hot  tenuous  mag¬ 
netospheric  plasma. 

There  are  also  ionospheric  signatures  of  the  LLBL  [Eastman  et 
aL,  1976].  These  signatures  are  controlled  to  a  great  extent  by 
the  plasma  flow  in  the  LLBL  which  is  mostly  in  the  antisunwatd 
direction,  although  a  region  of  relatively  slow  sunward  flow,  with 
occasional  large  velocity  peaks,  appears  to  be  present  towards  the 
magnetospheric  edge  of  the  layer  [Williams  et  aL,  1985;  Mitchell 
et  aL,  1987].  The  flow  is  mainly  perpendicular  to  the  magnetic 
field  and  therefore  creates  an  electric  field  across  the  LLBL  The 
resulting  electric  potential  differences  map,  at  least  approximately, 
along  the  magnetic  field  lines  to  the  ionosphere  where  they  drive 
horizontal  ionospheric  currents.  Where  the  ionospheric  electric 
field  changes  with  latitude,  the  above  currents  have  to  be  partially 
deflected  into  field-aligned  currents,  thus  forming  the  dayside  part 
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of  the  legion  I  current  system  [lijima  and  Potemra,  1976a,  b]. 
These  currents  flow  in  a  narrow  layer  into  (on  the  eveningside) 
and  out  of  (on  the  momingside)  the  LLBL,  where  they  are  again 
deflected  to  flow  across  the  magnetic  field  in  a  direction  opposite 
to  the  electric  field  there.  In  the  connecting  region  between  the 
ionosphere  and  the  LLBL,  field-aligned  potential  drops  are  likely 
to  occur.  In  the  current  system  described  above,  the  LLBL  acts  as 
a  magnetohydrodynamic  generator  while  the  ionosphere  plays  the 
roleofadissipator.  Recent  observations  [  Woch  etaL.  1993]  provide 
added  support  for  the  notion  that  the  LLBL  is  source  of  the  dayside 
region  1  currents. 

Compared  to  some  other  interfaces  in  the  solar  wind- 
magnetospheric  system,  relatively  little  theoretical  work  has  been 
devoted  to  the  LLBL.  Perhaps  the  first  qualitative  theoretical  model 
was  proposed  by  Coleman  (1 971].  Eastman  et  aL  (1976]  discussed 
the  LLBL  in  detail  with  emphasis  on  the  role  it  plays  in  transferring 
mass,  momentum,  and  energy  from  the  solar  wind  to  the  magneto¬ 
sphere.  In  that  work,  coupling  to  the  auroral  ionosphere  plays  an 
important  role  and  the  role  of  the  LLBL  as  a  generator  is  explained. 
Kan  and  Lee  (1980]  studied  imperfect  ionosphere-magnetosphere 
coupling  in  a  simple  evolving  boundary  layer,  using  concepts  of 
field-aligned  potential  drops  developed  by  Fridman  and  Lemaire 
(1980],  Chiu  and  Cornwall  (1980],  and  Lyons  (1980].  Their  model 
is  nonviscous.  the  result  being  that  the  ionospheric  drag  gradually 
slows  down  the  boundary  layer  plasma  as  it  moves  in  the  down¬ 
stream  direction.  The  role  of  viscosity  in  the  layer  was  studied  by 
Sonnerup  [1980]  in  a  one-dimensional  steady  state  model.  Lotko 
et  aL  [1987]  reexamined  this  model  by  including  a  simple  con¬ 
ductance  law  to  describe  the  relation  between  field-aligned  currents 
and  field-aligned  potential  drops.  In  those  papers,  induced  magnetic 
fields  were  not  included.  Wei  etal.  [1990]  examined  the  formation 
of  vortices  and  other  turbulent  structures  in  the  LLBL  caused  by 
the  Kelvin-Helmholtz  instability.  In  their  study  the  plasma  flow 
was  assumed  incompressible,  two  dimensional,  and  time  depen¬ 
dent  and  the  magnetic  field  was  taken  to  be  constant  Their  model 
included  viscosity.  Wei  and  Lee  [1993]  extended  that  model  to 
include  coupling  to  the  polar  ionosphere.  SiscoeetaL  [1991]  de¬ 
scribed  the  coupling  of  the  LLBL  model  of  Lotko  et  aL  [1987] 
to  a  high-latitude  boundary  layer  model,  and  Sitcoe  and  Maynard 
[1991  ]  also  included  coupling  to  the  region  2  current  system. 

Following  the  ideas  of  Sonnerup  [1980]  and  Lotko  et  aL 
[1987],  Phan  et  aL  [1989]  developed  a  fully  self-consistent  one- 
dimensional  model  of  a  viscous  LLBL  in  which  the  magnetic  field 
deformation  in  the  layer,  caused  by  the  currents  in  it,  is  included. 
The  present  model,  some  aspects  of  which  have  not  yet  been  nu¬ 
merically  implemented,  is  an  extension  of  the  analysis  by  Phan  et 
al.  to  allow  for  slow  variations  of  the  layer  in  the  flow  direction 
(the  negative  z  direction).  The  boundary  layer  occupies  the  region 
y  >  0.  starting  at  the  magnetopause  (y  =  0),  where  the  antisunward 
flow  velocity  has  a  maximum  value,  and  extending  earthward  into 
the  magnetosphere,  where  the  velocity  first  reverses  sign  and  then 
asymptotically  approaches  a  low  sunward  magnetospheric  value  as 
y  —  oo.  The  surfaces  *  =  ±Ho(x,y)  represent  the  northern 
and  southern  edges  of  the  layer,  where  the  magnetic  and  plasma 
pressures  reach  their  magnetospheric  values.  The  model  includes 
self-consistent  calculation  of  the  magnetic  field  in  the  layer  from 
the  currents  via  Ampfcre’s  law.  The  normal  component  of  the  cur¬ 
rent  is  continuous  at  the  surfaces  x  =  ±Ho,  whereas  the  jt  current 
switches  abruptly  to  zero  there;  for  |z|  >  Ho  the  current  is  field- 
aligned.  As  a  result  of  the  cross-field  currents,  the  field  lines  have 
approximately  parabolic  shapes  in  the  boundary  layer.  This  defor¬ 
mation  of  the  magnetic  field  along  with  the  deformation  caused  by 
field-aligned  currents  in  the  coupling  regions  between  z  =  ±Ho  and 
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the  ionosphere  are  important  features  of  the  model  because  they  pro¬ 
vide  the  mapping  along  LLBL  field  tines  from  the  equatorial  plane 
to  the  auroral  ionosphere.  However,  in  the  simplified  version  of  the 
model  to  be  discussed  here  the  lumped  properties  of  the  two  iono¬ 
spheres  and  the  low-pressure  (low  plasma  beta)  regions  between  the 
ionospheres  and  the  LLBL  proper  are  replaced  by  resistive  parallel 
plates  located  a  fixed  distance.  2Ho.  apart;  as  a  consequence,  the 
magnetic  mapping  is  not  included  in  a  self-consistent  manner  and 
field-aligned  potential  drops  are  excluded. 

The  paper  is  organized  as  follows;  The  mathematical  formulation 
of  the  model  is  presented  in  section  2  along  with  the  assumptions 
and  approximations  used.  Section  3  contains  a  description  of  the 
numerical  procedure.  In  section  4  we  present  three  runs  of  the  nu¬ 
merical  code,  and  in  section  3  we  summarize  the  results  and  discuss 
generalizations  of  the  present  model  that  remain  to  be  implemented. 

Development  of  Model  equations 

The  geometry  of  the  system  to  be  studied  is  shown  qualitatively 
in  Figure  I .  Figure  la  shows  the  dawns ide  LLBL,  viewed  from  the 
north  ecliptic  pole.  The  coordinate  system  is  as  follows;  the  y  axis 
is  across  the  LLBL  and  points  inward  toward  the  magnetosphere; 
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Fig.  1.  Boundary  layer  geometry.  (a)Dawnside  low-latitude  boundary  layer 
on  closed  field  lines,  viewed  from  above  the  north  ecliptic  pole.  Velocity 
profile  across  the  layer  is  shown  qualitatively,  including  a  legion  of  sunward 
flow.  Region  I  field-aligned  currents  flow  into  the  ionosphere,  (b)  View  of 
the  LLBL  from  the  Sun:  the  layer  extends  to  height  Ho  above  and  below 
the  equatorial  plane.  Parabolic  magnetic  field  lines  are  shown;  parabolas 
close  to  the  magnetopause  have  larger  curvature. 
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the  x  axis  is  directed  south  to  nonh  in  the  equatorial  plane;  and 
the  z  axis  is  perpendicular  to  the  previous  two  axes  and  points  in 
the  upstream  direction.  The  span  wise  width  of  the  LLBL,  which 
is  very  small  in  comparison  with  the  scale  of  the  curvature  of  the 
magnetopause  and  geomagnetic  held  lines  and  with  the  scales  for 
variation  in  the  stream  wise  direction  and  for  flaring  of  the  magne¬ 
topause  boundary,  permits  a  simple  mapping  of  the  slab  (z,  y,  a) 
coordinates  of  the  mathematical  model  described  below  onto  nat¬ 
ural  curvilinear  coordinates  such  as  (L,  M,  N)  boundary  normal 
coordinates.  In  Figure  1  a  velocity  profile  across  the  layer  is  shown 
qualitatively,  the  main  flow  being  in  the  -z  direction,  with  slow 
sunward  flow  in  the  magnetosphere.  Closed  magnetic  field  lines 
connect  the  LLBL  to  the  northern  and  southern  ionospheres.  In  the 
coupling  regions,  field-aligned  currents  occur.  Figure  lb  offers  a 
view  of  the  layer  from  the  Sun;  it  indicates  that  the  LLBL  extends 
to  a  height  Ho  above  and  below  the  equatorial  plane  in  the  x  di¬ 
rection  and  that  the  magnetic  field  lines  are  bent  into  approximate 
parabolas  with  increasingly  large  curvature  as  the  magnetopause  is 
approached.  In  the  current  version  of  the  model,  electrically  con¬ 
ducting  planes  at  *  m  ±Ho  with  appropriate  conductivity  are  used 
to  represent  the  ionospheric  substrate  and  the  low-pressure  regions 
between  it  and  the  effective  high-latitude  boundaries  of  the  LLBL. 

Boundary  Layer  Approximation 

The  plasma  flow  in  the  equatorial  LLBL  is  governed  by  the  steady 
state  mass  and  momentum  conservation  laws,  the  1  sen  tropic  law,  tire 
induction  law,  and  the  law  of  magnetic  flux  conservation; 
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nent  of  the  magnetic  field  is  of  the  order  of  (S/L)Bo  or  (6/H)Bq. 
As  mentioned  already,  another  consequence  of  the  boundary  layer 
approximation  is  that  the  magnetopause  surface  curvature  can  be 
ignored  if  the  characteristic  thickness  6  is  much  smaller  than  the  ra¬ 
dius  of  curvature,  ft,  a  condition  that  is  well  satisfied  for  the  LLBL. 
Thus,  neglecting  terms  of  order  6/ft  compared  to  terms  of  order 
unity,  one  is  allowed  to  use  a  Cartesian  coordinate  system  and  the 
flat  geometry  shown  in  Figure  lb.  By  neglecting  terms  of  order 
S2/L2,  62/HL .  and  S2/H2  compared  to  terms  of  order  unity  and 
by  setting  v,  =  0,  the  three  components  of  (2)  can  be  written  in  the 
following  approximate  form: 
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In  this  approximation,  inertia  and  viscous  forces  enter  only  the  z 
component  of  the  momentum  equation;  in  the  y  and  z  directions 
the  equations  express  a  static  balance  between  the  Vp  and  j  x  B 
forces.  The  motion  in  y  is  treated  purely  kinematically  via  the  mass 
conservation  law  (1).  The  viscous  term  in  (6)  has  been  assumed 
comparable  to  the  inertia  terms,  which  is  true  provided  the  viscous 
Reynolds  number.  Re  =  pV0L/ij,  is  of  the  order  of  L2/62;  over 
large  flow  distances,  one  expects  the  boundary  layer  thickness  to 
adjust  itself  to  satisfy  this  condition.  In  the  same  approximation, 
that  is.  again  neglecting  62/L2,  62/LH.  and  S2/H2,  the  currents  in 
the  LLBL  are  given  by  Ampferc’s  law; 


V  -  B  a  0  .  (5) 

Here  p,  p.  and  v  are  the  plasma  density,  pressure,  and  velocity, 
respectively.  Also,  r  is  the  viscous  stress  tensor.  B  and  mo  are  the 
magnetic  field  and  the  vacuum  permeability,  respectively,  and  7  = 
Cp/c.  is  the  ratio  of  specific  heats.  For  simplicity,  we  have  neglected 
resistivity  as  well  as  heat  conduction  and  viscous  dissipation  in  the 
LLBL.  These  effects  may  in  reality  be  of  some  importance  and 
should  ultimately  be  included. 

As  mentioned  already,  the  effects  of  the  two  ionospheres  are 
represented  by  parallel  resistive  plates  at  fixed  height  z  =  ±Ho- 
These  plates  are  assumed  impenetrable  in  the  present  version  of  the 
model  so  that  v,  =  0  at  z  =  ±Ho-  It  will  be  shown  in  the  next 
subsection  that,  in  the  approximate  equations  describing  the  LLBL, 
this  boundary  condition  implies  v,  =  0  in  the  entire  layer,  that  is, 
for  —Ho  <  z  <  Ho. 

The  boundary  layer  or  narrow-channel  approximation  is  now 
applied  to  the  above  equations,  taking  advantage  of  the  fact  that 
changes  across  a  narrow  layer  (in  the  y  direction)  occur  on  a  scale 
length  6,  which  is  much  smaller  than  both  the  scale  length  L,  for 
changes  along  the  main  flow  direction  (z),  and  the  scale  length  H . 
for  changes  in  the  third  direction  (z).  From  the  mass  conserva¬ 
tion  law  it  then  follows  that,  if  the  characteristic  velocity  in  the  z 
direction  is  Vo,  the  velocity  in  the  y  direction  is  of  the  order  of 
(S/L)Vo.  Similarly,  from  the  flux  conservation  law  it  follows  that 
if  the  characteristic  value  of  B.  and  of  ft*  is  Bo,  then  the  y  compo¬ 
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Series  Expansion  in  z 

We  proceed  now  to  expand  all  quantities  as  power  series  in  terms 
of  the  z  variable.  We  assume  a  boundary  layer  symmetric  with 
respect  to  the  equator  (the  plane  z  =  0)  and  express  each  dependent 
variable  in  terms  of  even  or  odd  powers  of  z,  according  to  whether 
they  are  symmetric  or  antisymmetric  with  respect  to  the  equatorial 
plane: 

Vj  =  »zo(*,|f)  +  v*:(z,y)—  +  ...  ; 

*>»  =  v*>(z<v)  +  vyi(z,v)-j^  +  ...  ; 
z  z3 

v*  =  v,x(z,y)—  +v,s  (*.*)jjjj +•••  (12) 


15 


2334 


Dkakou  rr  /u_  Low-Lxrrnjt*  Boundary  Urea  Mood. 


Bt  -  +  B*i(x,i)jp  +  ...  ; 

By  =  fl»i(*,y)^  +  B,s(*,y)-|j+  ■■  ; 

,2 

B.  =  B»o(*,y)  +  B,2(x,p)jp  +  ..  (13) 

,2 

*>  =  Po(*. y) +  *>:(*, t)-jp+  -  ; 

,2 

P  =  JH»(*,»)  +  Pl(x1|f)^j  +  ...  (14) 

In  the  model  we  shall  keep  only  the  lowest-order  terms  in  the  expan- 
sions  of  velocity  and  magnetic  field,  that  is.  terms  independent  of 
*  for  symmetric  quantities  and  terms  proportional  to  a  for  antisym¬ 
metric  quantities.  From  the  boundary  condition  v,  =  Oats  =  ±H0 
it  then  follows  that  «<(  =  0,  that  is.  that  o,  =  0  throughout  the 
boundary  layer.  We  emphasize  that,  in  contrast  to  the  nonevolving 
model  of  Phan  el  aL  { 1989),  the  lowest-older  terms  do  not  form  an 
exact  solution  in  the  present  model.  We  can  nevertheless  evaluate 
all  equations  at  *  =  0;  corresponding  quantities  are  denoted  by 
the  subscript  zero.  It  is  then  found  that  (3)  and  (8)  are  identically 
satisfied  and  that  (1).  (3),  (4),  (6),  and  (7)  become 

Jj(povio)  +  ^(Po^vo)  =  0  (15) 
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In  (16),  no  is  the  viscosity  coefficient  evaluated  in  the  equatorial 
plane,  and  in  (17)  the  quantity  Poo(x)  is  the  sum  of  the  plasma 
pressure  and  the  magnetic  pressure  in  the  equatorial  plane.  This 
total  pressure  is  independent  of  y  and  therefore  has  the  same  value 
in  the  boundary  layer  as  in  the  magnetosphere,  where  it  can  be 
considered  to  be  known  for  the  purposes  of  our  model.  By  use 
of  (15),  equation  (19),  which  represents  the  *  component  of  the 
induction  law,  can  also  be  written  in  the  form 

(20» 

The  x  and  y  components  of  the  induction  law  are  identically  satisfied 
at  x  =  0.  Finally,  note  from  (18)  and  (20)  that  po/pj  and  B^/po 
are  constant  along  streamlines. 

Equation  (8)  can  be  used  to  estimate  the  plasma  pressure  in  the 
LLBL  for  *  #  0,  after  the  magnetic  field  is  known.  This  equation 
implies  that  p  has  a  negative  quadratic  term  in  x,  in  addition  to  the 
term  po(z,  y).  which  describes  the  pressure  in  the  equatorial  plane. 
The  variation  in  pressure  described  by  this  term  is  consistent  with  the 


expectation  that  m  the  LLBL  the  plasma  pressure  should  fall  from 
a  high  value  in  the  equatorial  plane  toward  lower  values  at  the  high- 
latitude  edge  (x  =  ±£To)ofthe  LLBL,  although  the  magnetospheric 
pressure  level  is  not  actually  reached  in  the  present  model  due  to 
the  assumption  of  constant  Ho.  Extension  of  the  model  to  include  a 
variable  boundary  layer  height  so  that  the  plasma  pressure  reaches 
its  ambient  magnetospheric  level  at  x  -  ±Ho(x,  y)  is  planned,  as 
discussed  further  in  section  5.  It  is  also  noted  that  (8)  contains 
B this  field  component,  which  is  of  the  form  Bwi(x/H),  can  be 
calculated  from  Bt i  by  use  of  (5)  and  assuming  B, j(z,  j)  =  0 

In  general,  the  density  must  be  assumed  to  vary  with  x  as 
well,  falling  toward  a  lower  magnetospheric  value  at  z  =  ±H0. 
which  means  that,  in  addition  to  the  term  po(x,  y).  it  has  a  neg¬ 
ative  quadratic  term  in  x  as  well.  In  the  present  model,  where 
all  higher-order  terms  in  the  a  expansions  for  velocity  and  mag¬ 
netic  field  are  neglected,  the  quadratic  term  in  p  can  be  ob¬ 
tained  from  the  mass  conservation  law.  which  then  takes  the  form 
(vxod/dx  +  v»od/dy)(p2/0o)  =  0,  together  with  the  known  up¬ 
stream  density  distribution.  Thus  the  towest-otder  truncation  of 
the  x  expansion  of  velocity  and  magnetic  field  allows  us  to  calcu¬ 
late  second-order  terms  in  density  and  pressure  by  requiring  mass 
conservation  and  exact  force  balance  in  the  *  direction  to  be  main¬ 
tained.  If  the  expansions  of  v  and  B  are  carried  t*»  higher  order, 
the  lowest-otder  terms  in  these  quantities  and  the  terms  describing 
density  and  pressure  would  all  change.  Therefore  the  x  expansion 
used  is  not  exact,  in  the  sense  that  each  term  retained  is  not  the  same 
as  the  corresponding  term  in  the  expansion  of  the  exact  solution  in 
a  power  series  in  z;  rather  it  represents  a  polynomial  approximation 
to  the  exact  solution. 

The  system  of  (15).  (16).  (17).  (18).  and  (20)  contains  six  un¬ 
known  quantities,  namely  ®xo,  v^o.  B, i,  B, 0,  po.  and  po  Thus 
we  need  one  more  equation  to  form  a  closed  set.  This  equation 
will  be  provided  by  coupling  to  the  ionosphere,  as  discussed  in  the 
following  subsection;  it  will  lead  to  an  explicit  expression  for  Bt\ 
in  the  x  momentum  (equation  (16)). 

Ionospheric  Closure 

As  mentioned  already,  in  the  present  simplified  model,  the 
mapped,  lumped  properties  of  the  ionosphere  are  represented  by 
two  parallel  resistive  plates  at  fixed  height  x  =  ±Hq,  that  is,  at 
the  top  and  bottom  of  the  LLBL.  In  the  remainder  of  this  paper  we 
choose  the  characteristic  scale  height  Bin  (1 2M 1 4)  to  be  the  fixed 
boundary  layer  height  Ho.  The  net  current  flowing  on  the  surface 
of  the  plates  at  z  =  ±Ho  may  be  expressed  by 

JM  =  LEi  ,  Jy  =  ZEW  (21) 

where  J,.,  is  the  surface  current  density  and  £  is  an  effective 
height-integrated  conductivity  representing  the  lumped  response 
of  the  ionosphere  and  the  low-pressure  (low  plasma  beta)  regions 
between  x  =  ±H0  and  the  ionosphere.  The  components  of  the 
electric  field,  Ex  and  Ey,  imposed  on  the  plate  by  the  boundary 
layer  dynamo,  are  given  by 

Et  =  Et  =  -VyB,  \tm±Ht  ~  (6/l)E, 

Ey  —  Ey  |raiS|  —  VtB,  |«-±w0  (22) 

An  expression  for  the  conductivity  of  the  resistive  plates  that  repre¬ 
sent  the  ionosphere  in  our  model  was  given  by  Sormerup  [1980): 

Here  Ip  is  the  height-integrated  ionospheric  Pedersenconductivity. 
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Also,  £,  is  the  ionospheric  magnetic  field;  B  is  a  representative 
value  of  the  magnetic  field  *  component  in  the  equatorial  LLBL; 
dzjdz  is  the  ratio  of  a  length  element  dz,  in  the  ionosphere  to  a 
corresponding  length  element  dz  in  the  LLBL,  and  *  is  a  coupling 
factor  which  is  unity  when  the  magnetic  field  lines  that  connect 
the  LLBL  to  the  ionosphere  are  equipotentials.  In  its  present  form 
our  model  allows  for  field-aligned  potential  drops  in  the  coupling 
region  between  the  LLBL  and  the  ionosphere  only  in  the  average 
sense  obtained  by  letting  the  (actor  k  <  1.  This  restriction  will 
be  eliminated  eventually,  as  explained  in  section  5.  In  principle, 
the  effective  conductivity  £  is  a  (unction  of  r  and  y.  because  the 
mapping  factor  rfx./dx,  as  well  as  B  and  Ip  vary  in  space. 

We  now  impose  current  continuity  at  the  top  of  the  LLBL.  that 
is.  at  x  =  ±Ha,  which  leads  to 


.  .  dJt  dJu  dJu 

u  I-**-  +  "  if 


(23) 


In  this  expression.  J,  is  smaller  than  Jy  by  a  factor  of  the  order  of 
6/L\  thus  dJg/dz  is  smaller  than  dJt/d y  by  a  factor  62/L 2  and 
can  be  neglected.  Using  (1 1),  (21).  (22).  and  (23).  we  find 


Bz  |x«±//0  ~  -F#*o(£v»5«)  l*«±/fo  c(*) 


where  c(  x)  is  a  constant  of  integration,  and  the  negative  and  positive 
signs  on  the  right  side  of  the  equation  correspond  to  s  m  or 
z  =  -Hn.  respectively.  In  evaluating  the  left  side  of  the  above 
equation.  Bt j  and  all  higherterms  in  Bx  will  be  ignored,  as  before, 
so  that  Bz  =  Bz\(x/H0).  Similarly,  we  use  only  the  lowest-order 
terms.  Vzo  and  B, o.  in  evaluating  the  right-hand  side,  the  result 
being 


Bz i  =  -ncZvzoBzo  +  c(x)  (24) 


The  function  c(z)  is  given  by  the  boundary  conditions  at  y  =  oo 
(for  convenience,  the  interface  between  the  LLBL  and  the  magneto¬ 
sphere  is  thought  of  as  located  at  a  large  but  finite  y  value:  y  =  y^): 
at  that  location  t>io  assumes  the  low  magnetospheric  (sunward  or  an- 
tisunward)  value  vXOo(x).  the  magnetic  field  B, o  becomes  £«,(z). 
Bzt  becomes  Biiocfz).  the  conductivity  becomes  loo,  and  the  cur¬ 
rent  becomes  jy«>(z).  That  the  latter  quantity  can  be  nonzero 
indicates  that,  as  is  the  case  in  the  geomagnetic  tail  a  net  current 
may  be  flowing  across  the  magnetic  field,  from  the  LLBL  into  the 
magnetosphere,  or  vice  vena.  Equations  (10)  and  (24),  both  evalu¬ 
ated  at  y  =  oo.  can  now  be  combined  to  eliminate  Bzta>,  the  result 
being 


d  B 

c(z)  -  ftDl00v*coBoa-¥  ftoHojyao  +  Ho-^~  (25) 

Substitution  of  (24)  and  (25)  into  (16)  gives  the  following  final 
form  of  the  z  component  of  the  momentum  conservation  law  in  the 
equatorial  plane  (x  =  0) 


Po(VzO~  +  VyO^)Vzl,= 


+  (<Toot>«oB«  +i»»  +  --jf-)8* 

f).  (26) 

Note  that  <r,  as  defined  in  the  equation  following  (22),  can  in  prin¬ 
ciple  depend  on  x  and  y,  although  r  =  »„  =  const  is  used  in 
the  calculations  repotted  here.  Note  also  that  all  magnetospheric 
quantities  (denoted  by  subscript  oo)  are  functions  of  x.  in  general 
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Numerical  Procedure 

Finite  Difference  Scheme  and  Boundary  Conditions 

For  the  solution  of  the  five  coupled  equations  (15).  (17).  (18). 
(20).  and  (26),  of  which  four  are  partial  differential  equations  with 
independent  variables  r  and  y.  a  Crank-Nicolson  implicit  finite  dif¬ 
ference  scheme  was  implemented,  similar  to  marching  procedures 
that  have  been  used  in  the  past  for  fluid  mechanical  boundary  layers 
[e.g ..  Anderson  etai,  1984].  The  computational  domain  is  a  rectan¬ 
gular  region  in  the  xy  plane.  The  lines  y  =  Oandy  =  represent 
the  magnetopause  and  magnetospheric  boundaries  respectively.  At 
those  boundaries,  appropriate  boundary  conditions  are  imposed,  as 
described  below. 

At  y  s  0  the  velocity  tiio(i.O)  =  Vo(z)  is  specified  from  an 
appropriate  model  of  the  magnetosheath  flow  (located  in  the  region 
y  <  0).  assuming  Vzo  to  be  continuous  across  the  magnetopause. 
The  model  in  principle  allows  for  a  specified  plasma  flow  normal 
to  the  magnetopause,  but  in  the  results  repotted  here  we  have  used 
the  condition  vyo (x,  0)  =  0.  This  means  that  the  cold  plasma  mass 
flux  in  the  LLBL  remains  constant  as  one  moves  downstream  from 
the  computational  boundary  at  x  m  0  where  velocity,  magnetic 
field,  and  plasma  properties  are  specified  as  functions  of  y,  subject 
to  certain  consistency  requirements  to  be  explained  presently.  The 
restriction  Oyo(x,0)  =  0  also  means  that  the  upstream  station  can¬ 
not  be  located  at  local  noon  in  the  magnetosphere,  unless  a  delta 
function  mass  source  for  the  LLBL  is  assumed  to  be  present  there. 
In  principle,  the  initial  conditions  at  x  =  0  could  be  chosen  to 
match  data  from  a  satellite  traversal  of  the  LLBL;  the  evolution 
of  the  LLBL  downstream  from  the  satellite  crossing  could  then  be 
determined  by  the  model. 

At  y  s  yo«  the  x  velocity  Vzo.  the  magnetic  field  Bzo.  the  plasma 
density  and  pressure  are  fixed  to  their  specified  equatorial  magneto¬ 
spheric  values,  Vxoo(z).  Boa(z).  Poo(x)  and  foo(z).  respectively, 
and  are  assumed  to  approach  those  values  asymptotically  so  that 
(d/dy)tmWmm  2=  0.  For  reasons  of  mass  conservation  the  y  compo¬ 
nent  of  the  velocity  vyoo  cannot  be  specified:  it  is  determined  from 
the  numerical  calculations  and  represents  entrainment  of  magneto¬ 
spheric  plasma  (vyoo  <  0)  or  mass  feeding  of  the  magnetosphere 
from  the  LLBL  («,<»  >  0).  Finally,  at  the  downstream  boundary, 
located  at  x  =  L,  say.  no  conditions  are  needed  as  a  consequence 
of  the  parabolic  nature  of  the  problem. 

The  finite  difference  equations  comprise  a  non-linear  algebraic 
system  of  equations  that  involves  all  quantities  to  be  calculated 
across  n  grid  points  in  the  y  direction,  at  two  consecutive  steps 
in  the  -r  (main  flow)  direction.  As  mentioned  already,  initial 
conditions  for  various  quantities  are  imposed.  Note  that  po(0,  y )  and 
Bzo(0,  y)  must  be  chosen  to  obey  (17).  One  can  also  show  that  the 
initial  velocity  profile  Vyo(0,  y)  cannot  be  chosen  independently  but 
is  determined  from  the  other  variables.  After  B, o(0,  y),  po(0,  y), 
po(0,  y ),  and  Uxo(0,  y)  have  been  chosen.  tr*o(0,  y)  is  obtained  from 
the  equation: 

vxo(0,  y)  =  v„o(0, 0)  +  vzo  f  ~^r 
J 0  WzO 


[<rBio»*0-  (ffoofiooBoo  +  jyeo  +  — -'f”  }Bh> 


Mo  dz 


,  dP°°n  _  P0”k 

+  dz  K  _ 

IPO 


>vlo  v  d  ,  .1 


(27) 


Equation  (27)  is  derived  as  follows:  The  term  dvzo/dz  on  the  left 
in  the  x  momentum  equation  (equation  (26))  can  be  substituted 
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from  the  mass  conservation  law  (equation  ( 1 3)).  The  term  dpo/dz 
will  then  appear,  but  the  system  of  (17),  (18)  and  (20)  can  be 
solved  for  dpo/dz  as  well  as  dB, o/dz  and  dpo/dz.  By  further 
substitution  of  dpo/dz  into  (26),  only  derivatives  with  respect  to  y 
will  appear.  Equation  (27)  follows  by  first  forming  the  derivative 
(d/dp)(«yo/ vxo)  and  then  integrating.  It  is  evident  that  (27)  is  true, 
not  just  for  z  =  0.  but  for  any  value  of  z. 

After  the  problem  is  initiated  at  z  *  0  the  solution  is  found  si¬ 
multaneously  at  all  grid  points  across  the  boundary  layer  at  the  next 
step  in  the  main  flow  direction  (-*),  by  solving  the  nonlinear  alge¬ 
braic  system  mentioned  above;  this  solution  is  then  used  to  march 
forward.  The  system  is  solved  by  Newton’s  method  [Burden  and 
Fains,  1989],  which  is  a  fixed  point  (iterative)  procedure  requiring 
an  initial  guess  close  to  the  actual  solution.  Since  the  variations  of 
the  boundary  layer  in  the  z  direction  are  slow,  this  guess  is  provided 
by  the  values  at  the  previous  step  in  z. 

Asymptotic  Region 

It  is  one  of  the  goals  of  the  model  and  the  numerical  scheme 
to  be  able  to  predict  not  only  the  antisunwaid  LLBL  flow  but  also 
variations  of  slow  sunward  flow  observed  on  the  earthward  side  of 
the  LLBL.  It  should  be  noted  here  that,  partly  for  instrumental  rea¬ 
sons.  some  ambiguity  exists  in  the  data  regarding  the  latter  region: 
SckopkeetaL  [1981]  report  a  halo  region  adjacent  to  and  earthward 
of  the  LLBL.  where  the  density  and  temperature  are  intermediate 
between  those  in  the  LLBL  and  in  the  magnetosphere  and  where  the 
velocities  are.  for  the  most  part,  small  and  have  variable  direction. 
However,  sometimes  they  find  time  intervals  in  which  the  flow 
is  clearly  antiparallel  to  the  magnetosheath  flow.  Williams  et  aL 
[1983]  report  a  stagnation  region  adjacent  to  the  main  antisunward 
LLBL  flow  region,  where  velocities  are  small  and  have  variable 
direction,  and  a  region  of  small  but  steady  sunward  flow  earthward 
of  it  It  is  not  clear  whether  the  intermediate  region  contains  actual 
turbulent  flow  or  whether  it  is  simply  a  region  of  gradual  transition 
in  which  the  small  velocities  are  not  well-resolved.  As  described  in 
detail  below,  our  model  contains  an  asymptotic  flow  region  earth¬ 
ward  of  the  boundary  layer  proper  which  includes  a  flow  reversal 
followed  by  monotonically  increasing  and  then  constant  sunward 
flow  as  one  moves  from  the  boundary  layer  proper  into  the  magne¬ 
tosphere.  Inclusion  of  field-aligned  potential  drops  in  the  coupling 
region  between  the  LLBL  and  the  ionosphere  is  expected  to  change 
this  asymptotic  velocity  profile  into  one  that  includes  an  overshoot 
in  the  sunward  flow  adjacent  to  the  velocity  reversal  [Lotko  et  aL, 

mi). 

An  integration  procedure  dial  marches  forward  in  the  antisun¬ 
ward  direction  fails  within  a  few  steps  in  a  reverse  flow  region, 
a  difficulty  that  can  be  traced  to  the  term  »*o0c*o/dz  in  the  z 
momentum  equation.  This  term  retains  its  sign,  while  all  other 
terms  in  the  equation  containing  «*)  reverse  their  signs,  when  v*o 
reverses  sign.  In  the  reverse-flow  region  the  correct  marching  direc¬ 
tion  is  therefore  reversed.  This  problem  has  long  been  recognized 
in  computational  fluid  dynamics  and  different  solutions  have  been 
proposed.  Noticing  that  the  reverse  flow  is  usually  slow,  Reyhner 
and  Fliigge-Lotz  [1968]  suggested  neglecting  the  term  v^dv^/dx 
or  reversing  its  sign.  On  the  basis  of  that  idea  we  have  proceeded  to 
neglect  the  above  term  in  the  z  momentum  equation  in  the  reverse 
flow  region  as  well  as  in  a  narrow  region  on  the  magnetopause  side 
of  the  flow  reversal  point,  where  vro  is  negative  but  sufficiently 
small  so  that  vtodvto/dz  can  be  justifiably  neglected.  The  reason 
for  including  the  latter  region  is  that  there  is  also  a  second  difficulty 
at  the  flow  reversal:  when  the  reversal  point  falls  close  to  a  grid 
point,  a  numerical  singularity  appears.  We  emphasize  that  the  above 
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approximation  is  consistent  with  $  low-flow  convection  models  of 
the  inner  magnetosphere  (e  g.,  the  Rice  convection  model)  to  which 
the  LLBL  model  should  eventually  be  matched. 

To  implement  these  approximations,  we  have  terminated  our 
regular  computational  box  at  y  =  y»  but  appended  to  its  magne- 
tospheric  side  what  we  call  an  asymptotic  box,  y<»  >  y  >  y*.  in 
which  the  inertia  term  po«xodv«o/dz  is  neglected.  The  location 
y  =  jrt  is  different  in  every  step  of  the  calculation  along  the  -i 
direction:  in  the  tuns  presented  here  it  is  defined  by  the  require¬ 
ment  that  t>*o  at  that  location  be  9  km/s  or  less,  in  the  antisunward 
direction.  Furthermore,  in  the  asymptotic  box  we  take  pa.  pa.  and 
B>o  to  be  independent  of  y  so  that  pa  =  Poc(z).  po  =  Poo(z).  and 
B,o  =  Bao(x)  there.  With  the  above  assumptions  the  resulting  z 
momentum  equation  becomes 

dvia  .  .  „j  d* v*a  .... 

PoVfO~Q^~  —  —  ^oo(®*o  —  r>«)8go  +  too  •  (28) 

Equation  (28)  follows  from  (26)  by  setting  v*o dv^o/dx  equal  to 
zero  and  by  noticing  that  the  third  term  inside  the  parenthesis  of  the 
right-hand  side  of  (26)  cancels  the  magnetic  pressure  pan  of  the  term 
—dPoo(x)/dx.  The  plasma  pressure  pan  of  that  term  is  canceled 
by  jyeo  Boo  if  one  requires  the  solution  for  i>r0  to  be  independent  of 
y  as  y  —  oo.  Also,  the  viscosity  os  in  the  asymptotic  box  is  taken 
to  berjoo(x).  because  it  is  assumed  that  no  depends  in  a  general  way 
upon  po.  po.  and  B,a.  all  of  which  have  become  independent  of  y 
in  that  box.  If  the  y  variation  of  v^o  is  neglected  in  the  asymptotic 
box,  equation  (28)  has  a  simple  analytic  solution,  namely 

t»*o  =  vMoo(z)  +  [»*o(z,  y»)  -  Vxoo(z)]e_0<lK,'~*‘)  (29) 


-PoVpo{z,ya)+  v/pJvJoJz.y*) +  4700<t00S50 

“(*)= - iz -  (30) 


In  reality.  vpo  is  a  function  of  y  in  the  asymptotic  box;  the  above 
solution  is  therefore  only  an  approximation  to  the  actual  solution. 
In  the  runs  presented  in  this  paper,  the  above  solution  is  close  to  the 
exact  solution  when  v^q  has  a  slow  variation  with  y  as  it  usually 
does  away  from  the  upstream  location.  The  exact  solution  can  be 
found  by  integrating  equation  (28);  one  may  assume  that  Oyo  is 
approximately  known  from  the  previous  step,  so  that  (28)  becomes 
an  ordinary  differential  equation  for  t>*o- 

At  the  boundary  between  the  asymptotic  box  and  the  computa¬ 
tional  box.  y  =  y»,  we  require  that  v*o  and  v^>  be  continuous  and 
that 


«(*)v*o(*,y*)  + 


=  a(x)»*oo(x) 


which  guarantees  that  dv^o/dy  and  dvgo/dy  will  also  be  contin¬ 
uous  at  y  =  yt  Equation  (31)  is  obtained  from  (29)  by  differen¬ 
tiation.  As  can  be  seen  from  (29)  and  (31).  only  v^a  and  its  first 
derivative  can  be  made  continuous  at  y  =  y».  A  small  discontinuity 
in  the  second  derivative  of  v*o  remains  and  leads  to  a  corresponding 
discontinuity  in  the  first  derivative  of  the  current  j,.  In  the  results 
presented  in  this  paper  we  have  therefore  allowed  the  exponent 
o(z)  to  differ  slightly  from  t*e  value  given  by  (30),  in  such  a  way 
that  the  second  derivative  of  t>ro  is  also  continuous  at  y  =  y».  By 
differentiating  (29)  twice,  the  following  condition  is  found  at  that 
boundary: 


«J(z)t>*o(z,y»)  +  d  =  o:(x)v„„(x) 


Eliminating  o(z )  between  (3 1 )  and  (32)  yields  the  following  modi- 
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find  boundary  condition  at  the  interface  between  the  computational 
box  and  the  asymptotic  box: 


W  1  -•,•(*) 


(33) 


+  T-° 


'(,-S+'>S)"  = 


(39) 


We  can  finally  obtain  Vyo  in  the  asymptotic  box  by  integration  of 

the  mass  conservation  law: 


*«**,*)  -  '*>(«.») -(»- »>,”kj 


_dP~l*)  +  {RmVtoeBoo  +  Jtoo  + 


(40) 


1  d 

Poo(x)  rfx 

w) "  "  >)]  (34) 

Note  from  (34)  that,  in  general,  does  not  reach  a  constant  value 

but  retains  a  linear  variation  withy  as  y  —  oo. 

Since  the  nature  of  the  solution  in  the  asymptotic  box  it  such 
that  fjo  decays  at  a  relatively  slow  exponential  rate  to  e«(x),  the 
reduction  of  the  size  of  the  computational  box  that  results  from  use 
of  an  asymptotic  box  saves  a  great  deal  of  computer  time. 

Normalization 

We  now  explain  the  normalization  of  all  variables  in  (13),  (17), 
( 1 8),  (20).  and  (26).  which  leads  to  an  important  result  regarding  the 
role  of  viscosity  in  our  simplified  model.  The  normalized  variables, 
denoted  by  an  asterisk,  are 


The  viscous  Reynolds  number  has  been  scaled  out  of  the  equations 
because  of  our  particular  normalization  of  y  and  This  implies 
that  the  reference  value  of  the  viscosity  affects  our  system  only  as 
a  scaling  factor  for  the  width  of  the  LLBL  and  the  magnitude  of 
In  other  words,  a  thin  LLBL  with  low  viscosity  evolves  in  the 
same  way  as  a  proportionately  thicker  layer  with  higher  viscosity. 
However,  for  one  and  the  same  upstream  width,  boundary  layers 
with  different  viscosity  evolve  differently.  Note  also  that  we  have 
chosen  viscosity  to  depend  parametrically  on  pressure,  density  and 
magnetic  field,  as  indicated  by  the  factor  (p^'p^/flj J)  in  the  last 
term  of  (40).  For  example,  the  collisional  viscosity  transverse  to  a 
strong  magnetic  field  has  pt  =  -0.3.  pi  =  2.3,  pi  =  2  [Spitver, 
1962);  for  Bohm  diffusion  one  finds  pi  a  1 .  pi  =  0.  pj  =  1.  The 
reason  for  using  a  nonconstant  viscosity  is  that,  as  discussed  later 
on.  uniform  viscosity  throughout  the  LLBL  is  found  to  result  in 
excessive  entrainment  of  magnetospheric  plasma  and  high  accel¬ 
eration  of  this  plasma  downstream.  To  avoid  this  effect,  viscosity 
must  decrease  with  increasing  distance  from  the  magnetopause:  this 
is  accomplished, for  example,  for  positive  values  of  p(.  pi.  and  py. 


x*  =  x/ff*,  jr*  »  jn/fle/ffo. 


*  wv„,  *;  =  b; 


p*  =  po/pn,  P*  *  *00 

Pn*n 


iyoo  * 


(33) 


In  the  above,  ff0  is  the  characteristic  height  of  the  LLBL  while  V„ 
and  pw  are  reference  values  of  the  flow  velocity  and  density.  Abo, 
Re  it  the  viscous  Reynolds  number.  Re  =  pnVnH/i)n,  n*  beings 
reference  value  for  the  viscosity,  and  Rm  is  the  magnetic  Reynolds 
number,  defined  by  Rm  m  peoMl /mH*. 

The  dimensionless  system  of  (13),  (17),  (18),  (20).  and  (26)  now 
becomes  (the  asterisks  are  dropped  for  convenience): 

5^ (p**)  4  *  0  (36) 


p+^-ft.(*) 


(37) 


9  9 


o 


(38) 


Applications 

We  now  present  three  examples  of  LLBL  flow  obtained  from 
the  numerical  code.  We  first  describe  the  boundary  conditions  and 
initial  conditions  used  in  the  three  tuns;  in  die  remainder  of  the 
section  we  discuss  the  results  obtained. 

All  three  nuts  were  initialed  at  x  *  0.  The  initial  »*(»)  ve¬ 
locity  profile  is  a  somewhat  arbitrarily  chosen,  smooth  function  of 
p,  with  continuous  first  and  second  derivatives.  This  profile  sat¬ 
isfies  (31)  at  dx  interface  between  the  computational  box  and  the 
asymptotic  box;  in  the  latter,  it  includes  a  region  of  reverse  flow 
is  y  — >  goo  and  v,  becomes  independent  of  p  in  that  hmiL  With 
known  initial  profile  *,(p )  the  initial  v9  is  given  by  (27).  We  as¬ 
sumed  an  initial  number  density  decrease  of  80%  across  the  LLBL, 
from  10  protons/cm1  at  the  magnetopause  to  2  protons/cm3  in  the 
magnetosphere,  a  corresponding  increase  of  10%  in  the  x  compo¬ 
nent  of  the  magnetic  field,  from  40  to  44  nT.  and  a  constant  total 
pressure  of  1  nPa  across  the  layer.  The  initial  plasma  pressure  can 
then  be  calculated  from  (17),  the  result  being  that  the  temperature. 
7o(p)  S  po(p)/Apo(p).  where  R  =  ft/m,  is  the  gas  constant, 
varies  from  230  eV  at  the  magnetopause  to  773  eV  in  the  magneto¬ 
sphere.  The  parameters  chosen  for  these  runs  are  representative  of 
typically  observed  boundary  layer  crossings  [Sckopke  el  al,  1981; 
Eastman  and  Hones ,  1979].  They  result  in  an  average  electric  field 
component  across  the  layer  of  a  few  mV/m,  as  measured  directly  by 
Mover  [1984]  and  as  inferred  by  Mitchell  et  aL  [1987],  WUiams  et 
aL  [1985],  and  others  on  the  basis  of  particle  measurements.  Note 
that  the  magnetic  field  decrease  together  with  the  flow  acceleration 
lead  to  an  approximately  constant  value  of  this  electric  field  as  one 
moves  downstream  in  the  layer. 

The  values  of  the  characteristic  density  (p„  =  10  m,  kg/cm3) 
and  viscosity  i jn  ate  such  that  the  kinematic  viscosity  at  x  =  jr  = 
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0  is  t/n  s  i inlpn  =  10*  m2/s  [Sonnerup,  1980;  LaBelle  and 
Trtumann,  1 988]  in  the  fint  two  runs,  while  in  the  third  ran  viscosity 
has  been  reduced  by  a  factor  of  3  Also,  pi  . 2,  j  =  0  in  the  fint  ran  and 
pi  s  0.  pj,j  =  I  in  the  second  and  third  ran.  The  velocity  V»  used 
for  normalization  is  taken  to  be  280  km/s  and  Ho  =  10  Re  All 
calculations  ate  performed  with  the  ratio  of  specific  heats,  7  =  2. 
The  viscous  Reynolds  number  Re  used  in  the  fint  two  runs  is 
1.8  x  I04.  based  on  the  half  height  of  the  LLBL,  Ho.  and  on 
the  value  of  the  dynamic  viscosity  at  *  =  0,  p  =  0.  Note  that 
the  Reynolds  number  based  on  the  width  of  the  layer  is  at  least 
20  times  smaller.  The  magnetic  Reynolds  number  Rm,  based 
on  the  effective  conductivity  <Too  which  embodies  the  effects  of 
coupling  to  the  ionosphere,  is  0.1;  this  choice  leads  to  realistic 
values  for  the  field-aligned  currents  at  the  top  of  the  ionosphere, 
of  the  order  of  5  x  10“‘  A/m2.  To  achieve  that,  the  value  of  the 
effective  conductivity  is  reduced  by  a  factor  of  10  from  the  value 
5.96  x  10“*  mho/m  given  by  Sonnerup  [1980]  for  orange-segment 
mapping  and  perfect  coupling  to  the  ionosphere. 

The  boundary  conditions  are  as  follows; 

1.  At  the  magnetospheric  boundary,  that  is,  at  y  =  *«,,  vs  is 
constant,  vtoo  =  +10  km/s  in  the  sunward  direction,  while  the 
magnetic  field  falls  exponentially  in  the  downstream  (-z)  direc¬ 
tion.  according  to  the  formula  B,  =  £<»(*)  =  44exp(0.03c/ #£) 
nT.  This  simple  functional  form  is  not  necessarily  an  accurate  rep¬ 
resentation  of  the  actual  field  variation  in  the  magnetosphere,  but 
it  qualitatively  models  the  variation  of  Boo  with  downstream  dis¬ 
tance.  For  given  x  the  density,  plasma  pressure,  and  total  pressure 
are  constant  across  the  asymptotic  box  and  are  the  same  as  in  the 
magnetosphere.  They  can  be  calculated  from  the  following  three 
conditions: 

Bo.(»)  =  B«,(0) 

Poo(x)  Poo(0) 

Pe o(*)  _  P~(0) 

Pl>(*)  P«(9) 

p~(i)+^T  =  p~(l)  (4,) 

These  relations  follow  from  the  induction  law  and  the  isentropic  law, 
applied  along  the  dividing  streamline  passing  through  the  point 
v  =  yb  at  *  =  0,  and  from  the  y  component  of  the  momentum 
equation.  If  there  is  inflow  into  the  computational  box,  across 
y  —  n,  this  streamline  is  located  in  that  box;  therefore  the  values 
calculated  from  (41)  can  be  used  as  boundary  conditions  at  y  =  yi>. 
If  there  is  outflow  from  the  computational  box  across  y  =  jr»,  the 
dividing  streamline  is  located  in  the  asymptotic  box;  in  that  case 
the  values  from  (41)  can  be  used  as  boundary  conditions  at  y  =  ys 
only  if  plasma  leaving  the  computational  box  also  carries  with  it 
those  uniform  values  of  B,  p  and  p.  In  other  words,  a  region 
of  uniform  values  must  always  be  present  in  the  computational 
box  itself,  immediately  adjacent  to  y  =  g».  As  y  —  we 
require  a  balance  between  the  pressure  force  -dpa»/dx  and  the 
force  jyooBv,  from  which  the  current  jy«o(z)  that  flows  from  the 
LLBL  into  the  magnetosphere  can  be  calculated;  in  the  simulations 
presented  here  jyao  -A  0. 

2.  At  the  magnetopause,  that  is,  at  y  =  0,  the  plasma  is  as¬ 
sumed  to  accelerate  in  the  antisunward  direction  from  50  km/s,  at 
z  =  0,  to  290  km/s  at  z  =  -15  Re.  approximately,  150  km/s 
being  gained  within  the  first  5  Re.  as  described  by  the  formula 
wI(r,0)  =  -[308  -  258exp(x/5.75f?£)]  km/s.  Again,  this  for¬ 
mula  is  only  a  qualitative  representation  of  the  actual  variation  of 
v,(z,  0)  in  the  magnetosphere.  Since  we  do  not  allow  flow  across 
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the  magnetopause,  that  is,  since  t>„(z,0)  =  0,  the  line  y  =  0  is  a 
streamline.  Thus,  B,.  p.  and  p  at  y  =  0  are  defined  by  relations 
similar  to  (41). 

Figure  2  shows  results  of  the  first  simulation  example.  The 
bottom  row  of  the  figure  consists  of  three  plots  representing,  in 
order  from  left  to  right,  density,  temperature  and  *  component 
of  the  magnetic  field,  all  in  the  equatorial  plane,  as  functions  of 
y  at  three  different  locations  along  the  magnetopause,  namely,  at 
x  =  -0.25  Re.  -7.5  R  e  and  - 1 5  Re  The  pressure  profile  is  not 
shown  but  can  be  calculated  from  po  and  To:  it  has  the  property  that 
Po  +  Blo/2po  is  constant  across  the  layer.  The  z  and  y  components 
of  the  flow  velocity  and  the  z  component  of  the  current  at  the  top  of 
the  LLBL  are  plotted  in  the  middle  row  as  functions  of  y  at  the  same 
locations  downstream.  The  values  of  the  current  density  in  Figure 
2.  multiplied  by  B,/B, 0.  give  the  actual  field-aligned  currents  at 
the  top  of  the  ionosphere;  here  B,  =  5  x  10*  nT  is  the  ionospheric 
magnetic  field.  The  top  panels  consist  of  two  vector  plots  followed 
by  one  contour  plot  representing  the  z  component  of  the  current  The 
first  vector  plot  represents  the  actual  velocity  vectors;  in  the  second 
plot,  vy  was  magnified  by  a  factor  of  27  in  order  to  make  visible 
the  inflow  from  the  magnetospheric  boundary.  The  horizontal  axis 
in  all  three  top  panels  corresponds  to  the  y  direction  (across  the 
LLBL).  w  hereas  the  vertical  axis  is  the  distance  downstream  (the 
— z  direction).  The  width  of  all  panels  is  0.85  Re.  y  =  0  at  the  left 
being  the  magnetopause  boundary. 

The  dynamic  viscosity  in  this  first  example  is  taken  to  be  constant 
(i.e„  pi.2.3  =  0  in  (40))  throughout  the  LLBL  and  to  be  such  that  the 
kinematic  viscosity  at  the  magnetopause  and  at  z  =  0  is  equal  to 
the  reference  value.  10*  mJ/s.  The  profiles  in  Figure  2  indicate 
that  the  velocity  boundary  layer  is  becoming  thicker  and  that  its 
earthward  portions  accelerate  downstream,  as  plasma  enters  from 
the  magnetospheric  boundary  to  become  entrained  in  the  boundary 
layer  flow.  This  plasma  moves  into  the  LLBL  carrying  with  it  a 
higher  magnetic  field  and  lower  density.  Thus,  while  the  velocity 
layer  increases  in  width,  the  magnetic  field,  the  temperature  and 
density  layers  ail  decrease  their  widths,  as  shown  in  the  bottom  plots 
in  Figure  2.  Notice  also  the  decrease  in  density  and  temperature 
levels  with  increasing  z,  behavior  that  is  in  agreement  with  (41) 
and  with  our  choice  of  the  z  dependence  of  the  magnetic  field  in 
the  magnetosphere. 

The  plots  of  current  density  require  special  comment  The  t 
component  of  the  current  at  the  top  of  the  LLBL  is 

,  l  RBX 1  «./„  dvto  .  9B, o. 

+  <42) 

The  contour  plot  as  well  as  the  current  profile  show  that  the  current 
evolves  in  a  complicated  manner  creating  a  secondary  maximum  at 
some  distance  away  from  the  magnetopause:  the  two  terms  on  the 
right  in  (42)  have  opposite  signs  and  compete  in  a  way  that  results  in 
an  intervening  local  minimum.  As  seen  in  Figure  2,  this  effect  can 
even  create  a  current  reversal;  this  occurs  when  a  relatively  large 
negative  vt  makes  the  magnetic  field  profile  become  very  steep, 
thus  increasing  the  magnitude  of  the  negative  term  t?»o  dB^/dy  in 
(42).  We  draw  attention  to  this  point  because  the  observational  data 
often  show  a  highly  structured  behavior  of  the  region  1  currents  as  a 
function  of  latitude.  However,  field-aligned  potential  drops  M>n  in 
the  coupling  region  between  the  LLBL  and  the  ionosphere,  which 
were  included  by  Lotkoetai  [  1 987]  and  Phan  et  al  [1989).  artex- 
pected  to  modify  the  current  profile  significantly.  When  A4>n  =  0, 
the  value  of  the  current,  jt.  at  the  magnetopause  is  specified  by  the 
conductivity  £  and  by  the  magnetic-field  and  velocity  profiles,  as 
can  be  seen  from  (42):  in  the  runs  described  here  the  conductivity 
has  been  adjusted,  as  discussed  by  Sonnerup  (1980),  to  produce 
realistic  values  for  the  average  field-aligned  currents  at  the  iono- 
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Hg.  2.  Results  from  roa  1 :  coastaot  dynamic  viscosity.  Top  paadt  show  two  velocity  vector  plots  aad  a  coosow  plot  of  the  entreat 
density  j,  Ms  =  Ho-  Magnitude  of  v,  hat  bees  magnified  by  a  factor  27  in  the  eccood  vector  plot  to  show  the  e  atrainnirat  of 
magnetosphere  platina  into  the  layer.  Retnainieg  panels  thow  profile  plots  of  various  quantities  across  the  layer  at  z  =  -0.25  Rg, 
-7.9  Rg  and  -13  Rg  (narked  by  1,  2  and  3.  respectively).  The  computational  box  ends  at  y  =  0.18  Rg.  0.47  Rg.  and 
0.72  Rg  at  the  locations  1. 2  and  3.  respectively.  The  vM(y)  profiles  indicate  large  visctMteetmanmest  of  magaetotpheric  plasma 
and  a  rapidly  widening  velocity  layer.  The  current  density  profiles,  si  a  =  Hp  are  doable  peaked  with  sn  intervening  entreat 
reversal;.;*  has  an  absolute  max  imam  st  a  certain  a  value  (near  g  =  -7  Rg),  that  is,  at  a  certain  local  tune.  The  corresponding 
ionospheric  currents  are  given  by  the  currents  presented  here,  multiplied  by  a  factor  B,/Bbl  -  1667.  The  thickness  of  the  density, 
temperature,  end  magnetic  field  boundary  layers  decrease  at  the  plasma  moves  downstream,  even  thonghthe  thickness  of  the  velocity 
layer  increases. 


sphere  [tijima  and  Poiemra,  1976a].  On  the  other  hand,  when  first  increases  in  magnitude  with  increasing  flow  distance,  reach- 
A4>|l  ft  0  the  resulting  fourth  order  differential  equation  describing  ing  a  maximum  value  at  some  distance  |s|  downstream,  and  then 
the  LLBL  [Lotko  et  aL,  1987]  allows  one  to  specify  j,  =  0  at  the  reduces.  In  our  model,  two  effects  compete  to  produce  this  result, 
magnetopause  in  which  case  the  maximum  in  region  I  current  by  First,  the  flow  accelerates  downstream,  thus  increasing  the  velocity 
necessity  appears  at  some  distance  away  from  that  surface.  shear,  dv*o/<fy.  which  increases  the  magnitude  of  the  first  term  in 

The  longitudinal  variation  of  j,  shown  in  Figure  2  is  a  conse-  (42).  Second,  the  velocity  layer  increases  in  thickness  downstream 
quence  of  the  variation  of  the  flow  parameters  in  the  -z  direction,  as  a  result  of  viscosity,  thus  decreasing  the  shear,  and  consequently 
which,  unlike  the  previous  models  of  Lotko  et  aL  [1987]  and  Phan  the  magnitude  of  the  current  Initially,  the  first  effect  dominates  and 
et  aL  [1989],  continuously  modifies  the  two  terms  in  (42).  An  the  current  peak  increases  in  magnitude  as  one  moves  downstream, 
observed  feature  [lijima  and  Potemm,  1976a]  that  is  predicted  by  but  later  on  the  second  effect  takes  over  and  causes  the  maximum 
our  model  is  that  the  peak  intensity  in  the  region  I  current  density  current  to  decrease  again. 
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In  the  example  presented  above,  substantial  entrainment  of  mag¬ 
netospheric  plasma  into  the  LLBL  occurs.  As  shown  in  Figure  2.  the 
velocity  layer  becomes  thick  compared  to  the  density  layer  and  its 
eaithward  portions  accelerate  rapidly,  owing  to  the  assumption  that 
viscosity  remains  high  also  in  the  magnetospheric  plasma.  Obser¬ 
vations  by  Sckopke  et  aL  [1981]  and  Eastman  etaL  [1979]  do  not 
show  such  behavior  during  passage  from  the  magnetosphere  into 
the  LLBL  proper  the  velocity  increase  precedes  the  density  increase 
by  only  a  small  distance  so  that  only  a  small  amount  of  hot  mag¬ 
netospheric  plasma  is  entrained  by  the  tailward  flow  in  the  LLBL. 
Although  inclusion  of  field-aligned  potential  drops  in  the  coupling 
region  can  be  expected  to  modify  the  velocity  profile  predicted  by 
the  model  in  a  significant  way  [Lotko  etaL,  1987],  it  seems  plausi¬ 
ble  that  this  observed  weak  entrainment  of  magnetospheric  plasma 


is  accomplished  naturally  because  the  viscosity  varies  across  the 
layer.  In  the  next  simulation  example  the  viscosity  is  assumed  to 
drop  across  the  layer  in  proportion  to  the  drop  in  density  and  in 
inverse  proportion  to  the  increase  in  magnetic  field,  as  one  moves 
into  the  magnetospheric  plasma.  In  other  words,  pi  »  0.  pj  *  1. 
and  pj  =  I  in  (40). 

Figure  3  represents  our  second  example.  The  panel  arrange¬ 
ment  here  is  die  same  as  in  Figure  2.  Except  for  the  viscosity 
model  all  parameters,  and  the  initial  and  boundary  conditions  are 
the  same  as  in  the  first  example.  The  t»,  profiles  now  indicate  a 
curvature  reversal  at  the  approximate  location  of  maximum  slope  of 
the  density  and  magnetic  field  profiles,  a  result  of  the  variable  vis¬ 
cosity,  that  is.  of  the  fact  that  viscosity  drops  rapidly  as  one  leaves 
the  dense  boundary  layer  plasma  and  enters  the  tenuous  magneto- 


Fig.  3.  Results  from  run  2:  dynamic  viscosity  proportional  to  p/B,.  Format  is  the  same  as  in  Figure  2.  The  computational  box 
ends  at  v  —  0.18  Ag,  0.22  Ae.wdO.33  A e  at  the  locations  1, 2.  and  3,  respectively.  The  velocity  boundary  layer  is  now  much 
thinner.  Notice  also  the  curvature  reversal,  near  the  magnetopause,  of  the  profiles  of  v*  (y)  at  the  location  of  maximum  slope  of  the 
profiles  of  p(y)  and  B  ,(y ).  The  vv(y)  profiles  show  evidence  of  the  partial  relaxation  of  die  system  from  its  initial  conditions.  The 
secondary  current  density  maximum  is  now  higher  than  the  one  at  the  magnetopause  and  current  densities  are  higher  overall. 
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spheric  plasma.  The  field-aligned  current  again  develops  a  reversal 
region  followed  by  a  secondary  maximum  at  some  distance  away 
from  the  magnetopause,  which  is  now  higher  than  the  value  at  the 
magnetopause. 

The  first  u,  profile,  at  *  =  -0.23  Re,  differs  in  shape  from 
vv  profiles  further  downstream.  This  behavior  is  caused  by  a  rapid 
relaxation  from  the  arbitrary  initial  velocity  profile  vz(0,  y ):  within 
a  short  distance  downstream  the  flow  settles  down  to  a  more  regular, 
slow  evolution  determined  mainly  by  the  boundary  conditions.  As 
an  illustration  of  the  fading  memory  of  the  system  with  respect  to 
the  initial  velocity  distribution,  we  refer  to  the  analytic  solution 
given  by  Sonne rup  [1980].  With  constant  values  of  Vo.  Brf,  and 
<t,  and  with  dPao/dx  =  0,  vIOO  =  0  and  >voo  =  0,  (26)  has  the 
asymptotic  solution 


Vx0=  -Voe-""",  vy0  =  0  (43) 


as  *  — ►  —oo,  where  S„  is  the  viscous  length,  used  by  Atkinson 
[1967],  Sonnerup  [1980],  and  others: 


6U  =  ( 


loo  .«/» 
OooBl*’ 


(44) 


This  solution  also  satisfies  (15),  (17/,  (18),  and  (20)  identically. 
The  density  profile  can  be  arbitrary  in  this  solution.  Any  chosen 
initial  v*(0,y)  profile  and  its  associated  u„(0,  y)  profile  evolves 
in  such  a  way  that,  after  some  relaxation  distance  along  the  main 
flow  direction.  vx(z,  y)  becomes  identical  to  (43)  and  vv(x,y)  = 
0.  In  this  case  the  system  loses  memory  of  the  initial  velocity 
state  completely,  the  final  state  being  determined  entirely  by  the 
boundary  conditions.  However,  information  about  the  initial  density 
profile  remains.  In  the  general  case,  memory  of  nonconstant  initial 
magnetic  field,  density,  and  temperature  profiles  is  retained  and  a 
state  that  is  independent  of  z  is  never  reached,  unless  the  boundary 
conditions  at  y  —  0  and  y  =  y<»  become  independent  of  z. 

In  the  first  example  the  relatively  higher  viscosity  at  the  location 
of  the  secondary  maximum,  reduces  the  velocity  shear,  thereby 
spreading  the  velocity  boundary  layer  and  decreasing  the  magnitude 
of  the  first  term  in  (42).  The  result  is  relatively  lower  magnitudes  of 
the  current  at  that  location  and  overall.  In  the  second  example  the 
velocity  boundary  layer  is  narrower,  and  the  shear  is  larger  at  the 
location  of  the  secondary  maximum,  as  a  result  of  lower  viscosity  in 
the  magnetosphere.  The  secondary  current  maximum  is  larger  and 
lower  |t>y|  values  are  generated,  that  is,  less  magnetospheric  plasma 
is  being  entrained  into  the  layer.  The  total  potential  difference  across 
the  LLBL,  between  the  magnetopause  and  the  location  where  the 
velocity  reverses  sign  in  this  run,  is  1.65  kV  initially,  4.4  kV  at 
a  distance  IS  Re  downstream,  and  5.65  kV  at  a  distance  15  Re 
downstream.  These  values  of  the  potential  drop  across  the  LLBL  are 
consistent  with  typically  measured  potential  drops  [Mozer,  1984], 

In  our  last  simulation,  shown  in  Figure  4,  we  reduce  the  reference 
viscosity  by  a  factor  of  3,  from  vn  =  10*  to  vn  =  3.3  x  10*  mJ/s 
but  maintain  its  dependence  on  p  and  B  (pi.j  =  1,  jn  =0). 
As  a  result,  the  viscous  Reynolds  number  increases  by  a  factor 
of  3.  We  initiate  this  case  with  the  same  upstream  conditions  as 
in  example  2  and  use  the  same  boundary  conditions  as  well.  At 
the  downstream  locations,  z  =  -7.5  Re  and  r  =  — 15  Re,  the 
velocity  boundary  layer  is  now  thinner  than  in  the  previous  two  runs. 
There  is  a  relatively  large  initial  outflow  of  plasma  (v„  >  0)  across 
the  magnetospheric  boundary  during  which  the  initial  excessive 
amount  of  plasma  participating  in  the  downstream  motion  leaves 
the  layer.  As  a  result  of  lower  viscosity  close  to  the  magnetopause, 
no  reversal  of  field-aligned  current  occurs. 

In  Figure  5  we  plot  magnetic  field  lines  from  run  2  at  z  = 
-7.5  Re  in  three  zz  planes,  namely,  at  the  magnetopause  plane. 


y  =  0.  and  two  planes  parallel  to  it,  located  at  y  =  700  km  and 
y  =  4700  km.  respectively.  The  view  is  horn  outside  the  dawnside 
LLBL,  looking  towards  the  magnetosphere  and  with  the  Sun  on  the 
right.  Figure  5  shows  that  the  field  lines  are  bent  into  approximate 
parabolas;  their  curvature  is  greatest  near  the  magnetopause  and  is 
least  at  the  magnetospheric  boundary;  some  curvature  remains  at 
y  =  4700  km  because  jtv>  yt  0;  as  a  result  of  the  increasing  flow 
speed  and  decreasing  magnetic  field  the  curvature  increases  as  one 
moves  downstream.  Notice  that  the  scale  in  the  z  direction  has 
been  exaggerated:  at  z  =  10  Re  the  tines  are  in  reality  displaced 
~  1.5  fis  sunward  of  their  intersection  point  with  the  equatorial 
plane.  We  believe  that  in  the  real  LLBL  such  a  field  configuration 
may  reduce  or  eliminate  turbulence  in  the  layer,  because  magnetic 
flux  tubes  of  different  curvature  do  not  easily  interchange  their 
locations.  The  magnetic  field  configuration  described  above  and  its 
stabilizing  effect  were  also  noted  by  Southwood  [1979]. 

Discussion 

In  this  paper  we  have  presented  initial  results  of  a  steady  sate 
numerical  model  of  the  low-latitude  boundary  layer  on  closed  field 
lines.  The  model  is  experimental:  our  purpose  is  to  see  if  it  can 
predict  certain  features,  such  as  spatial  distribution  of  field-aligned 
currents  Rowing  into  the  ionosphere,  that  are  observed  in  the  data 
and  to  examine  the  effect  of  different  parameters  in  the  model  such 
as  the  viscosity  formula,  viscous  and  magnetic  Reynolds  numbers, 
upstream  conditions  and  boundary  conditions  at  the  magnetopause 
and  in  the  magnetosphere.  In  contrast  to  the  models  by  Phan  et 
at.  [1989],  Lotko  et  aL  [1987],  and  Sonnerup  [1980]  the  present 
model  includes  variation  of  boundary  layer  properties  such  as  flow 
velocity,  magnetic  field,  plasma  pressure,  and  density  as  one  moves 
in  the  downstream  (-z)  direction,  thus  allowing  for  evolution  of 
the  LLBL  in  the  flow  direction.  The  variation  of  all  quantities  with 
distance  (z)  away  from  the  equatorial  plane  is  parameterized  in  a 
simple  way.  The  flow  evolution  in  the  z  direction  is  governed  by 
a  balance  between  inertia  forces,  j  x  B  forces,  pressure  forces, 
and  viscous  forces,  whereas  in  the  two  perpendicuiardirections  (y 
and  z)  the  boundary  layer  approximation  results  in  satic  balance  of 
forces.  The  currents  in  the  LLBL  are  calculated  in  a  self-consistent 
manner,  via  Ampfere’s  law.  As  a  result,  the  magnetic  field  com¬ 
ponents  Bi  and  B.  are  allowed  to  be  of  comparable  magnitude, 
the  field  lines  being  approximately  parabolic,  with  vertices  point¬ 
ing  in  the  antisunwatd  direction  and  with  maximum  curvature  for 
field  lines  adjacent  to  the  magnetopause.  At  the  magnetospheric 
edge  of  the  layer  the  curvature  is  smaller  but  need  not  be  zero  as  a 
consequence  of  allowing  a  finite  cross-field  current  jy«e(x)  to  flow 
from  the  boundary  layer  into  the  magnetosphere.  We  distinguish 
here  between  the  Phan  et  aL  [1989]  model,  where  the  magnetic 
field  lines  are  exact  parabolas,  and  the  present  model  where,  owing 
to  the  (slow)  z  variation,  they  are  only  approximate  parabolas. 

We  emphasize  the  following  results  obtained  from  the  present 
version  of  the  model:  (1)  The  velocity  boundary  layer  increases 
in  thickness  downstream  as  a  result  of  viscous  entrainment  leading 
to  inflow  into  the  boundary  layer  region  of  low  density  magneto- 
spheric  plasma  carrying  higher  magnetic  field.  This  inflow  leads 
to  thinning  and  steepening  of  the  density,  temperature  and  mag¬ 
netic  field  profiles.  (2)  The  velocity  profile  tends  to  relax  from 
an  arbitrarily  assumed  upstream  shape  to  shapes  that  are  mainly 
governed  by  the  boundary  conditions.  However,  the  LLBL  never 
loses  memory  of  the  initial  magnetic  field,  density  and  pressure 
profiles.  (3)  The  thickness  of  the  velocity  profile  is  greatly  influ¬ 
enced  by  viscosity  which  enters  the  system  as  a  scaling  factor  (see 
equation  (35)),  larger  viscosity  resulting  in  thicker  velocity  profiles. 
(4)  The  field-aligned  currents  that  flow  into  the  ionosphere  from  the 
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Fig.  4.  Rentas  from  ran  3:  dynntac  viacoatay  propoitiooal  to  »/B,  and  three  tiraa  lower thn  in  ran  2.  The  companion*!  box 
ends  at  g  «  O.H  Ag.  0.1 3  A*,  and  0.21  Ag  at  the  locations  1,2  and  3,  respectively.  The  velocity  layer  is  tinnaer  dun  in  runs  I 
•ad  2.  NocaneaercvernlO*  <  0)  is  preaeat  io  this  case. 


upper  and  lower  edges  of  the  LLBL,  at  *  =  ±H0,  represent  the 
dayside  region  1  current  system  [lipma  and  Potemm,  1976a],  or 
at  least  those  portions  of  it  that  we  generated  on  closed  field  lines. 
These  currents  are  found  to  form  a  secondary  maximum  —ay  from 
the  magnetopause.  Between  the  magnetopause  maximum  and  the 
secondary  maximum  a  current  minimum  occurs,  often  including  a 
reversal  of  the  current.  The  region  1  field-aligned  current  peak  is 
found  to  reach  a  maximum  at  a  certain  local  time,  as  observed  by 
lijima  and  Potemra  (1976a).  However,  experience  with  the  lotto 
et  aL  (1987]  model  indicates  that  significant  changes  in  the  details 
of  velocity  and  field-aligned  current  distributions  win  occur  when 
field-aligned  potential  drops  are  incorporated. 

In  the  future  the  model  will  be  improved  in  the  following  two 


major  respects.  First,  the  ionosphere  will  be  incorporated  in  a  more 
realistic  way.  The  connection  region  between  the  LLBL  and  the 
ionosphere  is  a  narrow  channel  starting  at  the  upper  and  lower  edges 
of  the  boundary  layer,  that  is,  at  *  =  ±Ho,  where  the  currents 
are  field-aligned  and  extending  along  the  magnetic  field  into  the 
ionosphere.  In  this  region  the  magnetic  field  will  be  calculated 
self-consistently,  rather  than  by  using  the  constant  mapping  factor 
(dx,/dx)  employed  here.  Also,  field-aligned  potential  drops  will 
be  included  as  in  the  works  by  Lotto  et  aL  [1987]  and  Phan  et  aL 
[ 1 989],  according  to  the  empirical  formula  j||,  =  A'(d«~d*)-  Here 
>11,  is  the  field-aligned  current  into  the  ionosphere,  A  is  an  effective 
field-aligned  conductance  density,  and  and  #,  are  the  electric 
potential  at  x  =  ±Ho  and  in  the  ionosphere,  respectively.  The 
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7.5  -x/Re  6 


Fig.  5.  Magnetic-field  lines  in  the  diwuide  LLBL  from  ran  2  at  *  = 
-7.5  Re.  View  from  outside  the  layer,  looking  towards  the  magnetosphere, 
with  the  Sun  on  the  right  The  lines  have  approximately  parabolic  shape 
(exaggerated  in  the  plot  by  a  factor  of  3.3  in  the  horizontal  scale)  and 
maximum  curvanne  near  the  magnetopause. 


ionosphere  will  continue  to  be  treated  as  a  conductive  substrate; 
combining  Ohm’s  law  with  current  continuity  and  the  boundary 
layer  approximation,  we  can  show  that 


The  height-integrated  Pedersen  conductivity  may  have  a  spatial 
dependence  or  a  dependence  on  the  electron  precipitation  associated 
with  the  field-aligned  potential  drops. 

The  second  improvement  of  the  model  is  to  include  a  variable 
height  H0{z ,  y)  of  the  LLBL,  to  be  calculated  self-consistently  with 
the  other  variables.  In  this  case  the  velocity  component  v,  will  not 
be  zero,  except  at  z  =  0;  it  will  be  allowed  to  have  a  linear  a 
dependence.  The  surface  z  =  ±Ho(x,  y)  will  be  loaded  where 
the  plasma  pressure  in  the  layer,  as  determined  by  the  z  momentum 
equation  (equation  (8)).  has  dropped  to  its  magnetospheric  value. 
By  requiring  the  velocity  to  be  tangential  to  that  surface  it  can  be 
shown  that  the  mass  conservation  law,  at  z  =  0,  becomes 

9  9 

— (ifopofxo)  +  q^(HoPov*o)  —  0 

As  a  further  generalization,  this  equation  may  be  modified  to  allow 
plasma  entry  from  the  upper  and  lower  edges  of  the  LLBL,  at 
*  =  ±£fo(*,»)- 

Further  improvements  of  the  model  can  be  made  by  permitting 
mass  diffusion  across  the  magnetic  field. 
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APPENDIX  2:  EXTENSIONS  OF  THE  MODEL 


In  this  chapter,  the  equations  are  presented  that  describe  the  LLBL  with 
(a)  variable  height,  2Ho(x,y),  in  the  z  direction  and  (6)  field-aligned  potential  drops 
in  the  force-free  coupling  region  between  the  LLBL  and  the  southern  and  northern 
ionospheres. 

2.1  Variable  Boundary  Layer  Height 

The  basic  formulation  of  the  low-latitude  boundary  layer  model  remains  as 
described  in  Appendix  1,  but  in  addition,  the  half-height,  H0  is  now  a  function  of  x 
and  y  to  be  calculated  self-consistently.  It  is  assumed  that,  at  the  location  z  =  ±H0, 
the  plasma  pressure  has  fallen  to  its  corresponding  magnetospheric  value,  poo(i,  z), 
which  is  considered  to  be  a  known  function  of  x  and  z  only.  At  the  same  location,  a 
force-free  (jxB  —  0)  region  begins  that  connects  the  LLBL  to  the  two  ionospheres. 
Equations  (l)-(5)  describe  the  plasma  flow  in  the  equatorial  LLBL.  By  using  the 
series  expansions  (12)  and  (14),  the  mass  conservation  law,  equation  (1),  evaluated 
at  z  =  0,  takes  the  following  form 
d(po«*o)  ,  d(p0vv0)  p0vxl 

<45) 

Notice  that  vx  is  now  non-zero  everywhere  in  the  LLBL  except  at  z  =  0.  It  is  required 
that  the  plasma  velocity  vector  be  tangential  to  the  upper  and  lower  boundary  surfaces 
at  z  =  ±Ho". 

{17  •  V(z  -  tfo(z,y))}*=±//o  =  0 


or,  by  keeping  the  zeroth  and  first-order  terms  only  in  the  series  expansion  in  z, 


VxodHo  VyodHo  vz\Hq  ~ 


dx 


H 


=  0 


(46) 
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d(poHoVgo)  t  d(poHovvo)  _  n 
dx  +  dy 


By  substituting  vMi/H  from  (46)  into  (45),  the  mass  conservation  law  finally  becomes 

(47) 

The  induction  law  should  now  be  used  in  the  form  of  equation  (19),  instead  of 
(20),  since  the  latter  was  derived  by  use  of  (15),  which  is  no  longer  valid. 

The  x  and  y  components  of  the  momentum  conservation  law,  equation  (2), 
evaluated  at  z  =  0,  are  given  by  equations  (16)  and  (17).  In  the  boundary  layer 
approximation,  described  in  Appendix  1,  the  2-momentum  equation  is  given  by 

(48) 


dz 

In  the  above  equation,  the  magnetic  terms  on  the  right  are  of  order  {H/L)  compared 
to  the  term  on  the  left-hand  side,  and  the  inertia  and  viscous  terms  are  of  order 
(H/L)3,  where  H  and  L  are  the  characteristic  scale  lengths  for  changes  in  the  z 
direction  and  the  x  direction,  respectively.  At  this  point,  one  can  proceed  to  different 
levels  of  approximation:  if  H/L  is  of  order  unity,  one  would  need  to  solve  the  full 
2-momentum  equation;  if  H/L  is  much  less  than  unity,  one  can  either  neglect  the 
terms  of  order  (H/L)2  and  keep  the  terms  of  order  H/L,  or  neglect  both.  In  the 
following  we  will  neglect  terms  of  order  (H/L)3,  a  choice  that  is  consistent  with  the 
assumption  that  second  order  terms  in  the  z  expansions  could  also  be  neglected. 
According  to  this  latter  assumption,  terms  of  order  H/L  should  be  kept;  however, 
significant  simplifications  result  from  neglecting  those  terms  as  well.  The  two  different 
options  for  the  2-momentum  equation  are  stated  below. 

(1)  Terms  of  order  H/L  and  (H/L)2  are  neglected 

In  this  case,  the  2-momentum  equation  (8)  reduces  to  the  simple  pressure  balance 

Bl  ,  , 

P+2^  =  Po(*,y)  (49) 
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where  po(x,  y)  is  the  plasma  pressure  in  the  equatorial  plane.  If  the  left-hand  side  of 
the  expression  is  now  evaluated  at  z  =  ±f/0  where  the  pressure  has  become  equal  to 
the  magnetospheric  plasma  pressure,  Poo(x,  //<»),  evaluated  at  the  top  of  the  boundary 
layer,  one  finds 

=  0  (50) 

According  to  the  series  expansion  in  z,  the  magnetospheric  plasma  pressure  contains 
a  term  proportional  to  z3,  i.e., 

z  ^ 

Poo(x,  z)  =  pOQ(x,  0)  -  (jj)  pa(x) 

At  z  =  Hoo(x)  this  pressure  has  reached  the  value  that  is  assumed  to  be  present 
at  the  top  of  the  LLBL,  z  =  Ho,  throughout  the  layer,  as  shown  in  Figure  6.  Its 
value,  Poo(x,/f<x>)  is  considered  known  so  that  equation  (50)  can  be  used  to  calculate 
H0(x,y). 


(2)  Terms  of  order  (H/L)  retained 


In  this  case,  the  z-momentum  equation  (8)  reduces  to  the  static  force  balance 


"( 


Bx i  dBzo  Byi  dBx o  \  Hq 


(51) 


2p0H 3  \  po  dx  '  p0  dy  )  2  H 

which  replaces  (50)  as  the  equation  for  Ho{x,y).  In  this  expression  the  y  component 
of  the  magnetic  field  can  be  obtained  from  equation  (5),  which  simplifies  to  the  form 


dB: 


*x  ,  dBvi  _ 


=  0 


(52) 


dx  dy 

provided  dBz/dz  is  zero  which  is  the  case  if  quadratic  and  higher-order  terms  in  the 
series  expansion  of  Bx  are  neglected. 

The  system  of  equations  (46),  (47),  (16),  (17),  (51),  (18),  (19)  and  (52)  contains 
nine  unknown  quantities,  namely,  v*o,  Vyo,  u,!,  BxU  ByX,  Bx 0,  p0,  po  and  H0.  If  option 
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(1)  above  is  used  By\  does  not  appear,  except  in  equation  (52),  and  equation  (50) 
replaces  (51).  To  obtain  a  closed  system  of  equations,  coupling  to  the  ionosphere 
must  be  incorporated.  We  first  apply  current  continuity  at  z  =  ±H0-  The  surface 
z  =  H0(x,y)  is  defined  by  the  equation  F(x,y,z)  =  z  -  H0(x,y)  =  0;  the  unit  vector 
normal  to  this  surface  is  n  =  V In  what  follows,  all  quantities  just  below 
the  surface  F  =  0  carry  the  superscript  —  while  quantities  just  above  F  —  0  carry 
the  superscript  +.  Current  continuity  across  the  surface  z  —  H0(x,y)  implies  that 

i"*n  =  i+*n  (53) 


where  the  current  j~  at  z  —  Hq  is  given  by  (9)-(ll).  By  use  of  these  equations  and 

the  magnetic  field  as  given  by  (13)  one  then  finds 

1  dBMdHo  1  dH0  Btl  dB*  ,  H0  dBxl  +dH0  ,  ,+  dH0  .+ 

/io  dy  dx  +  no  dyK  H  dx  *  n0H  dy  Jx  dx  +  3*  dy  3x  {  } 

The  force  free  condition  above  the  surface  F  —  0,  namely  j+  x  B+  =  0,  along  with 

B~  =  B+,  which  implies  that  there  are  no  surface  currents  at  z  =  Ho,  result  in 


_ 


HB 


*o 


a  _  hbm 


(55) 


it  HqBx\  '  -  H0Byl 

If  the  magnetic  field  lines  in  the  coupling  region  are  equipotentials,  the  relation¬ 
ship  between  the  electric  field  Ey  at  the  top  of  the  LLBL  and  2?,  in  the  ionosphere  is, 
to  lowest  order,  given  by 


dy_ 

'dyi 


Ei  =  Ey-Zr*vt0Bzo 


dy_ 

dyi 


(56) 


Here,  dyi  is  the  ionospheric  counterpart  of  a  length  element  dy  in  the  equatorial 
LLBL;  a  relation  between  these  lengths  is  obtained  from  the  conservation  of  flux  in 
a  magnetic  flux  tube 
Bm  dx 


dyi  _ _ _ 

dy  Bi  dii 


(57) 
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where  B,  is  the  ionospheric  magnetic  field  and  dx,  is  the  ionospheric  counterpart  of 
an  equatorial  length  element  dx.  Note  that  if  the  equatorial  vector  elements  dx  and 
dy  are  orthogonal,  the  corresponding  ionospheric  elements  dx,  and  dyx  are  generally 
not  orthogonal  and  vice  versa.  However,  in  the  boundary  layer  approximation  this 
feature  of  the  mapping  does  not  enter  the  equations  explicitly.  The  ratios  dx/dxi  and 
dy/dy ,  are  known  only  after  the  field-line  geometry  in  the  coupling  region  has  been 
calculated  self-consistently,  as  explained  later  in  this  chapter.  In  the  simplest  case, 
one  may  assume  a  known  average  value  for  this  ratio,  as  was  done  in  Appendix  1. 

In  the  ionosphere,  the  height  integrated  current  is 

I,  =  E p(Ei  +  vnx  Bi )  -  E h(E{  +  vn  x  B{)  x  z,  (58) 


In  the  above  equation,  E /  is  the  height-integrated  Pedersen  conductivity,  Ej/  is  the 
height  integrated  Hall  conductivity,  and  t>„  is  the  ionospheric  neutral  wind  velocity. 
By  applying  current  continuity  at  the  top  of  the  ionosphere  and  by  assuming  constant 
vn  and  B,  we  can  use  the  boundary  layer  approximation  to  find 

4  =  ™ £<**>  =  ££<&•*>  (59) 

where  y'u.  is  the  field-aligned  current  leaving  the  ionosphere.  From  charge  and  mag¬ 
netic  flux  conservation  it  also  follows  that 


4  _  B< 

3t  B«> 

By  use  of  (56)  and  (57)  equation  (59)  becomes 


Bx  dii  d  ( 

4  “  R  _  Am  /).,  ^ 


E  p  ~r~vroBzo 


) 


Bm  dx  dy  \  '  Bz o  dx 
Finally,  by  use  of  (55)  and  (60)  equation  (54)  gives 


1  dBrtdHo  1  dH0  Bxi 

Ho  dy  dx  no  dy  H 


dBz o  Ho  dBxi 
dx  }  fioH  dy 


(60) 


(61) 
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(62) 


ill,  ( HoBs i  dHo  t  H0B vi  dHo  D  ^ 

W~JT~d^  +  ~H~W  ~ 

In  the  approximation  (2)  in  which  terms  of  order  H/L  are  retained,  all  terms  in  (62) 

must  be  used;  in  the  approximation  (1)  in  which  both  terms  of  order  H/L  and  ( H/L )2 

are  neglected  compared  to  terms  of  order  unity,  it  simplifies  to 

1  dBypdHo  1  dH0  Btl  dB*  H0  8Bxl  _  .  B* 

Ho  dy  dx  Ho  &V  H  &x  HoH  dy  ^ 

In  this  case,  equations  (61)  and  (63)  close  the  system  (46),  (47),  (16),  (17),  (50),  (18), 

(19)  and  (52)  mentioned  above,  with  one  more  dependent  variable,  namely  j||(.  For 

convenience,  the  complete  system  of  equations  to  lowest  order  is  repeated  here: 

d(p0H0vx0)  d(poH0vvo)  _  n 
dx  +  dy 

,  d  ,  d  ,  dP^x)  ,  1  D  D  ,  9  dvx 

Mv-a9i +  — ir~ +  + 

d  d  .  po  . 

UlOo - H  VyO -TT-)—  =  0 

dx  dy  pi 

d(B mVxq)  d(BMavvo)  _ 
dx  +  dy 

VxodHp  _  VypdHo  VxiHp  _ 

dx  dy  H 

dBxl  dBvi  __ 
dx  +  dy  " 


2.2  Field-Aligned  Potential  Drops 


In  addition  to  all  features  of  the  LLBL,  the  ionospheric  substrate,  and  the  force- 
free  region  connecting  the  two,  described  in  the  previous  sections,  the  field-aligned 
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current  in  the  ionosphere,  j\\ . ,  is  now  assumed  to  be  given  by  the  lumped  relation 

j\U  =  *(&  -  <f>i)  (64) 

where  k  is  am  effective  field-aligned  conductance  density  and  the  subscripts  e  and 
t  denote  quantities  evaluated  at  the  top  of  the  boundary  layer  and  the  ionosphere, 
respectively.  The  potential  distribution  in  the  LLBL  region  is 


<{>'  =  -(  Evdy  =  -  /  v^Btody 
Jo  Jo 


(65) 


The  potential  distribution  in  the  ionosphere  is  related  to  the  ionospheric  electric  field 
by  E{  =  —  By  use  of  this  expression  for  Ei,  and  by  noting  that  in  the  boundary 
layer  approximation  EXi  <  Eyi,  equation  (59)  becomes 


.  _  dy  9  tv  dy  dd>i\ 
^  dy,  dy  Pdy{  dy 


(66) 


By  substituting  <j>i  from  (64)  into  (66)  and  by  using  (57)  and  (65),  one  then  finds 


(67) 


Equation  (67)  instead  of  (61)  now  completes  the  system  of  equations  derived 
above.  Note  that,  without  field-aligned  potential  drops,  i.e.,  with  it  =  oo,  equation 
(67)  reduces  to  equation  (61). 


2.3  Asymptotic  Solution 

The  system  of  equations  in  the  previous  sections  cannot  be  integrated  in  the 
region  of  slow  sunward  flow  by  a  marching  procedure  in  the  —  x  direction,  for  the 
reasons  mentioned  in  Appendix  1.  As  in  the  approach  taken  in  that  Appendix,  the 
term  v^dv^/ dx  in  the  x-momentum  equation  will  be  neglected  in  the  slow  sunward 
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flow  region.  However,  it  is  not  possible  to  derive  a  simple  solution  there  under  the 
same  assumptions  as  in  Appendix  1,  namely,  that  £?*> ,  p  and  p  are  independent  of 
y.  The  reason  is  that  the  2-momentum  equation,  in  its  simplest  form  (i.e.,  equation 
(50)),  forces  the  height  H0  to  be  a  function  of  both  x  and  y.  The  remaining  equations 
are  then  inconsistent  with  those  assumptions  and  without  them,  we  have  not  found 
a  simple  straightforward  method  to  reduce  the  equations  describing  the  LLBL  in 
the  asymptotic  region  to  a  system  of  ordinary  differential  equations.  One  then  has 
the  option,  either  to  integrate  this  system,  without  the  above  inertia  term  in  the  x- 
momentum  equation,  using  the  same  numerical  method  as  in  the  computational  box 
in  Appendix  1,  or  to  make  the  assumption  that  the  height  of  the  layer  is  a  function 
of  x  only  and  does  not  vary  with  y.  In  the  latter  case,  the  z-momentum  equation 
cannot  be  used  to  satisfy  the  condition  that  the  pressure  drops  to  its  magnetospheric 
value  at  z  =  H0.  Instead,  the  pressure  above  and  below  the  surface  z  =  H0  will 
be  different,  except  in  the  magnetosphere  where  Hq  =  H^x)  is  chosen  to  provide 
continuity  of  pressure  across  this  surface.  As  explained  in  the  next  two  sections,  one 
may  account  for  this  pressure  difference  by  including  a  surface  current  at  z  =  Ho ,  as 
shown  in  Figure  6.  Note  that  this  pressure  difference  has  been  ignored  in  the  earlier 
models  by  Sonnerup  [1980]  and  Lotko  et  al.  [1987]. 

(1)  No  Surface  Currents  at  z  =  Ho 

In  this  approach  the  resulting  pressure  at  z  =  Ho  differs  from  the  magnetospheric 
pressure;  the  difference  can  be  considered  as  the  pressure  exerted  mechanically  by 
the  surface  itself  on  the  LLBL;  it  results  in  an  inconsistent  evaluation  of  the  currents 
and  the  magnetic  field  above  the  surface  z  =  Ho-  This  approach,  however,  leads 
to  a  simple  analytic  solution  in  the  asymptotic  region,  under  certain  assumptions  to 


be  explained  presently.  It  should  be  mentioned  that  the  analysis  presented  below  is 
essentially  the  same  as  that  performed  by  Lotko  et  al.  [1987]  and  used  by  Siscoe  et 
al.  [1991],  and  Siscoe  and  Maynard  [1991]. 

As  in  Appendix  1,  we  will  assume  that  Bt 0,  p  and  p  are  independent  of  y.  This 
assumption  is  not  necessary,  but  it  is  now  allowed  by  the  equations,  and  it  is  used 
to  illustrate  the  simplest  case.  With  the  above  assumptions,  the  current  continuity, 
equation  (63),  becomes 

H0  dBxi  Bzo  .  /co. 


Po  H  dy 


D.  -7|l i 


and  the  x-momentvm  equation  becomes 

dPoo(x)  1  D  _  dvxo 

*-§? - dx  +  =  ’ 

Equation  (67)  can  be  written  in  the  following  form  if  one  assumes  that  E p,  By  and 
dxi/dx  are  also  independent  of  y  in  the  asymptotic  region. 

(70) 

where,  using  the  notation  of  Lotko  et  al.  [1987] 

0=|l^,  A  =JIE  (71) 

Bz o  dx  V  « 

We  now  operate  on  (69)  with  (1/c2  —  X2d2/dy2)d/dy  and  express  dBx\jdy  in  terms 
of  j||.  from  (68).  Furthermore,  we  neglect  the  ^-dependence  of  Vyo  in  the  asymptotic 
region.  The  result  is 

,i  x2  d2  x  d  (  d2vx<>  ,  ,dVx0\  b]q  i  x2d2 n 

(?  - A  a^Ty  r  ■ w~ _  = 0  (72) 

which,  after  substitution  of  the  right  most  term  from  (70)  and  one  integration,  be¬ 


comes 


d_  (  d2 _ 1_\  (  d_  _  poVyojx^X 

dy  \dy2  c2X2)  \dy  r)Q  ) 


,  E  pBIq  _  E  PB% 
Ul0+  rioX'iHoB-0*0  rtoX'HoB?* 
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The  above  equation  has  the  general  solution 

Vgo  =  u*oo  +  e~ctv  (Acos(c2y)  +  Bsin(c2y))  (74) 

where  the  growing  exponential  has  been  discarded.  The  constants  A  and  B  are  to 
be  determined  from  the  boundary  conditions  at  y  =  y&.  Furthermore,  c\  and  c2  can 
be  obtained  by  requiring  the  first  and  second  derivative  of  to  be  continuous  at  the 
interface  between  the  computational  and  the  asymptotic  region. 


(2)  Surface  Current  at  z  =  Ho 


In  addition  to  the  assumptions  used  in  case  (1),  a  surface  current  K  is  now  assumed 
to  flow  in  the  surface  z  =  Ho,.  This  current  exerts  a  vertical  force  on  the  plasma 
in  the  boundary  layer  which  represents  the  net  effect  of  jy  currents  flowing  in  the 
“triangular”  (shaded)  region  in  Figure  6,  and  accounts  for  the  difference  in  pressure 
at  the  edge  of  the  LLBL  (at  z  =  Ho  =  #<»(*))  and  the  magnetospheric  pressure. 
This  requirement  is  given  by  the  relation 


Ky - - - 


(75) 


where  the  superscripts  -j-  and  —  denote  quantities  evaluated  above  and  below  the 
surface  z  =  Ho,  respectively.  In  the  above  equation,  the  magnetic  field  in  the  surface 
z  =  Hoo  is  the  average  between  the  value  above  and  below  this  surface.  In  the 
boundary  layer  approximation,  the  boundary  condition  on  the  magnetic  field  there 
requires  that  1?+  —  B~  =  jioKv,  where  B~  =  (Ho/H)Bxl .  The  plasma  pressure  below 
the  surface,  p~,  which  can  be  taken  from  the  ^-momentum  equation  (49),  is  equal 
to  po  —  BllHo/2poH2]  the  plasma  pressure  above  the  surface  is  the  magnetospheric 
plasma  pressure  at  z  =  Hoo,  i.e.,  it  is  Poo(x,Hoo)-  With  those  substitutions,  the 
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relation  (75)  gives  the  surface  current  Kw 


^^2  D]  i  t%PO  Poo 

■'+  ~7r~ 


(76) 


The  condition  of  current  continuity  at  z  =  Ho  should  now  include  this  surface  current, 
i.e.. 


.  •+  *  ,  dKy 

3  -n  =  y*n  +  -~ 
ay 


(77) 


which  leads  to 


Ho  dBxi  Bx o .  dKy 

=  ~Wn- +  W 


(78) 


Equations  (69),  (70)  and  (78)  contain  x  only  as  a  parameter;  they  are  in  effect  ordinary 
differential  equations  which,  along  with  (76),  can  be  used  to  calculate  v*o,  Bxl  and 
j||.  as  functions  of  y  at  any  x  location. 


2.4  Self-Consistent  Coupling  Region 

It  is  an  important  objective  of  a  boundary  layer  model  to  provide  the  actual 
mapping  of  the  magnetic  field  lines  from  the  equatorial  LLBL  to  the  ionosphere.  In 
order  to  accurately  provide  this  information,  the  model  must  include  a  self-consistent 
calculation  of  the  magnetic  field  deformation  in  the  coupling  region,  i.e.,  it  must 
incorporate  the  mapping  factor  dx/dii  in  a  self-consistent  manner.  This  can  be 
accomplished  by  an  iterative  procedure,  as  follows:  In  the  first  iteration  a  constant 
value  of  the  above  ratio  is  used,  equal  to  the  value  given  by  an  internal  magnetospheric 
model,  such  as  the  Tsyganenko  model.  After  the  model  equations  have  been  solved 
numerically,  the  magnetic  field  geometry  in  the  LLBL  will  be  known.  A  force-free 
boundary  layer  model  can  then  be  used  to  calculate  the  further  deformation  of  the 
magnetic  field  in  the  connecting  region  that  is  caused  by  field-aligned  currents  there. 
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This  calculation  will  provide  the  new  mapping  factor  dx/dii  to  be  used  in  a  second 
iteration  in  the  LLBL  model,  and  so  on.  A  self-consistent  boundary-layer  model 
of  the  connecting  region  is  currently  under  development  by  Professor  Lotko  and  his 
students  [see  Appendix  5]. 
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APPENDIX  3:  NUMERICAL  METHOD 


In  Appendix  1,  it  was  explained  how  the  system  of  equations  (15),  (17),  (18), 
(20)  and  (26),  i.e.,  the  mass  conservation  law,  the  y  component  of  the  momentum 
equation,  the  isentropic  law,  the  induction  law,  and  the  x  component  of  the  momen¬ 
tum  equation,  respectively,  all  evaluated  at  the  equatorial  plane,  (z  =  0),  were  solved 
by  a  numerical  procedure.  The  reader  is  reminded  that  the  model  represented  by 
those  equations  does  not  include  field-aligned  potential  drops  in  the  coupling  region 
between  the  LLBL  and  the  ionosphere,  and  assumes  that  the  LLBL  has  constant 
height  in  the  z  direction. 

As  already  mentioned,  the  computational  box  is  a  rectangular  region  in  the  xy 
plane.  The  lines  y  =  0  and  y  =  yoo  represent  the  magnetopause  and  the  magneto- 
spheric  boundaries,  respectively;  appropriate  boundary  conditions  are  applied  there. 
At  the  line  x  =  0,  the  upstream  conditions  across  the  boundary  layer  (i.e.,  in  the  y 
direction)  are  imposed.  The  box  is  divided  into  two  parts:  the  first  part,  attached 
to  the  magnetopause,  extends  from  y  =  0  to  y  =  yy,  in  this  part,  the  velocity  in  the 
negative  x  direction  is  higher  than  a  chosen  small  positive  value,  vmtn.  In  the  present 
chapter,  the  finite  difference  method  used  in  this  part  of  the  box  is  explained  in  some 
detail.  The  analytical  method  used  in  the  box  attached  to  it,  called  the  asymptotic 
box,  and  extending  from  y  =  yb  to  y  =  yoo,  was  explained  in  Appendix  1.  A  schematic 
of  the  computational  region  and  the  asymptotic  region  is  shown  in  Figure  7. 

The  system  of  equations  contains  five  unknown  quantities,  namely,  v*o(x,y), 
vyo{x,y),  Bzo(x,y),  p0{x,y)  and  Po(x,y);  we  wish  to  find  the  x  dependence  (along 
the  main  antisunward  flow  direction)  and  y  dependence  (across  the  layer)  of  these 
variables.  The  system  is  written  in  a  finite  difference  form  according  to  the  Crank- 
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Nicolson  scheme,  to  be  explained  presently.  There  are  n+2  grid  points  in  the  y 
direction,  including  the  two  boundary  grids,  and  the  calculation  can  proceed  for  as 
long  as  needed  along  the  x  direction.  The  equations  are  expressed  in  backward  differ¬ 
ence  form  in  the  x  direction;  this  allows  all  variables  to  be  calculated  (simultanuously 
at  all  n  grid  points  in  the  y  direction)  at  only  one  step  in  the  x  direction  each  time, 
based  on  the  values  of  the  previous  step.  The  grids  are  equally  spaced  in  the  x  and  y 
directions,  but  Ax  ^  Ay,  and  also,  care  is  taken  so  that  Ax/ Ay  «  v^/v^.  The  latter 
condition  guarantees  that  information  about  the  velocities  in  the  x  and  y  direction 
travel  a  distance  of  about  one  step  in  the  corresponding  direction  in  time  At,  i.e., 
A t  =  Ax/vxo  «  Ay/vyo-  This  marching  procedure  in  the  downstream  direction  takes 
advantage  of  the  slow  parabolic  evolution  of  the  flow;  in  fluid  mechanics  it  has  been 
found  to  be  the  most  economical  method  as  far  as  computational  time  is  concerned. 

At  every  step  in  the  x  direction,  the  number  of  grid  points,  n,  along  the  y 
direction  is  different  in  general.  Before  the  calculation  at  each  step  in  x  starts,  a  simple 
check  on  the  magnitude  of  v*o  at  the  right  boundary  is  performed;  for  computational 
convenience  this  check  is  applied  at  the  previous  step.  If  the  value  is  lower  than 
vmin ,  as  many  grid  points  as  necessary  are  discarded  at  the  right  edge  of  the  box  and 
the  calculation  at  the  present  step  is  performed  on  the  remaining  grid  points.  If  the 
value  is  higher  than  umi„,  grid  points  are  added:  the  values  of  all  variables  at  those 
new  points  at  the  previous  x-step  are  taken  from  the  asymptotic  (analytic)  solution 
explained  in  Appendix  1.  It  is  sufficient  to  check  the  value  of  v^o  &t  the  very  last 
grid  point  on  the  magnetospheric  edge  of  the  numerical  box  because  the  velocity  is 
expected  to  change  with  y  relatively  slowly  at  that  boundary. 
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3.1  Finite  Difference  Method 


The  finite  difference  form  of  the  x-momentum  equation  in  non-dimensional  form 
is  given  below.  The  superscript  i  denotes  the  grid-point  number  in  the  x  direction 
while  the  subscript  j  denotes  the  grid-point  number  in  the  y  direction.  In  the  following, 
the  notation  for  v*o,  v^o,  Bto,  Po,  jy  and  po  is  also  changed,  for  convenience,  to  u,  t>, 
B,  p,  j ,  and  p,  respectively.  The  x-momentum  equation  becomes 

!«(»;+■<') +(i- »)(^«s)i(uf*  -  -i) , 

A* 

-  <4ti) + (i  - 

2Ay 


eUnU^+j^  +  ^j  +  (1  „Bco  +  >,»  +  Bj+ 

+  Rn(l  -  *)(uB)j+ 


{*  ((‘ft'./ASl  -  “)*')  -  -  “/-•>)  + 

(1  -  *)  (/‘j+i/alxj+i  -  “J)  -  .  *  —  0, . . . i  =  1. —  >n  (79) 


In  the  above,  the  viscosity  is  assumed  to  be  of  the  form  p  =  pl,1pp2/Bp3.  Also,  the 
notation  Pj+1/2  =  (Pj+i  +  Pj)/ 2  and  Pj-1/2  =  (Pj  +  Pj- 1  )/2  is  used,  with  0  <  9  <  1, 
where  9  is  a  weighting  factor.  For  9  =  0  the  method  is  explicit;  in  this  case,  the  von 
Neumann  stability  constraint  presents  a  severe  limitation  on  the  marching  step  size 
[e.g.,  Anderson  et  al.,  1984].  For  9  =  1  the  method  is  fully  implicit  with  truncation 
error  0( Ax)  -I-  O(Ay)2.  The  value  for  9  used  here  is  1/2,  which  gives  the  Crank- 
Nicolson  implicit  scheme  with  truncation  error  0(Ax)2+O(Ay)2,  when  all  coefficients 
and  properties  are  evaluated  at  the  expansion  point  (i+1/2,  j).  In  this  case,  no 
stability  constraint  arises  from  the  von  Neumann  analysis. 
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The  y-momentum  equation  is  not  a  differential  equation.  One  could  possibly 
eliminate  one  of  the  variables  and  not  use  this  equation.  However,  for  the  purpose 
of  keeping  the  non-linear  terms  as  simple  as  possible  in  the  code,  the  equation  was 
written  in  a  discrete  form  and  was  included  in  the  system  along  with  the  other  four 
differential  equations.  In  non-dimensional  form,  this  equation  is  written  as: 

p)+'  +  =  (ft.)*1  (80) 


The  mass  conservation  law  is  expanded  at  the  point  (i+1/2,  j-1/2)  as  follows: 

Pj  ui  Pjuj  +  Pj-iuj-i  Pj-iuj-t 
Ax 

,  Pi  VJ  Pj-lvj-l  +  Pjvj  ~  Pj-lVj-l  _  A  ,ai 


Notice  that  the  derivative  across  the  layer  (in  the  y  direction)  is  written  in  a  back¬ 
ward  difference  form  here,  while  the  derivatives  across  the  layer  in  the  x-momentum 


equation  are  in  a  central  difference  form.  As  a  result,  no  boundary  condition  on  vy  is 
required  at  the  right  edge  of  the  computational  box. 

The  isentropic  law  and  the  induction  law  are  written  in  a  form  that  contains 


convective  terms  similar  to  the  inertia  terms  in  the  x-momentum  equation;  therefore, 
they  are  expressed  in  a  similar  finite  difference  form.  It  should  be  noted,  however, 
that,  in  b^th  equations,  the  derivative  across  the  layer  is  not  represented  exactly  by  a 
central  difference,  as  in  the  x-momentum  equation,  but  as  a  combination  of  backward 
and  forward  differences,  with  weighting  factors  0\  and  (1  —  $i),  respectively.  Without 
those  weighting  factors,  a  numerical  oscillation,  which  grows  as  the  computation 
proceeds  downstream,  appears  between  alternate  grid  poins  in  the  values  of  B,  p  and 
p.  Empirically,  the  value  6\  =  0.35  was  found  to  work  well. 
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The  isentropic  law  and  the  induction  law  have  the  following  finite  difference 


form. 

[guf '  +  (l  -  0)u;l  ((p/ !>•’)•+'  -  (,/pX 

Ax 

.  „  M-’r  -  &/?)&) + (i  - «.)  (wpX\  -  wpX') 

+  '  2Ay 

.  „  m  (<P/P’)}  -  (p/p’)}-,)  +  (1  -  «>)  ((P/P’)}+.  -  (P/P’) j) 
+(1  -  - — - 


(82) 


For  a  given  i  value,  equations  (79)-(83)  comprise  a  system  of  5xn  algebraic 
equations  with  5xn  unknowns,  namely,  u^+1,  vj+1,  Bj+1,  />}+1,  and  pj+1,  where  j  = 
l,...,n.  Boundary  values  at  the  magnetopause,  i.e,  at  j=0,  are  required  by  the 


numerical  procedure  for  all  five  quantities.  At  the  grid  point  n+1,  i.e.,  at  the  boundary 
between  the  numerical  and  the  asymptotic  box,  boundary  values  are  also  required  for 
B ,  p  and  p.  Boundary  value  for  v  is  not  required.  Note  that  B,  p  and  p  at  both  the 


magnetospheric  boundary  and  the  magnetopause  boundary  should  be  consistent  with 
the  relations  (41)  and  relations  similar  to  (41),  respectively,  as  explained  in  chapter 
2.  The  boundary  condition  on  u  at  the  grid  point  n+1  is  a  mixed  condition,  i.e., 
it  contains  both  the  value  and  the  two  first  derivatives  of  the  quantity,  see  equation 
(33).  This  equation  is  implemented  numerically  as  follows  [Smith,  1984]:  equation 
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(84) 


(33)  is  written  at  the  grid  point  n+ 1,  i.e.,  at  y  =  as 
un+2  -  2«n+i  +««=  for  all  i 

«*+l  “  «oo 

Equation  (79)  is  then  written  at  the  grid  point  n+1,  and  the  value  of  un+a  is  substi¬ 
tuted  from  (84).  This  results  in  one  additional  equation  with  one  additional  unknown 
quantity,  namely  un+*.  Therefore,  the  final  algebraic  system  consists  of  /Xn)+1 
equations  with  (5xn)+l  unknowns.  The  method  used  to  solve  this  non-linear  system 
is  described  in  the  next  section. 


3.2  Newton’s  Method 

Newton’s  (or  Newton-Raphson’s)  method  is  one  of  the  most  powerful  and  well- 
known  numerical  methods  for  solving  a  root-finding  problem  /(x)  =  0  [ Burden  and 
Faires ,  1989].  To  indroduce  Newton’s  method,  assume  that  the  function  /  is  twice 
continuously  differentiable,  and  xa  is  an  approximation  to  the  root,  r,  of  /(x)  =  0 
such  that  /'(x0)  ^  0,  where  f  is  the  derivative  of  /.  One  may  then  consider  the 
Taylor  expansion  of  /(x)  around  the  point  x„: 

/(x)  =  /(x0)  +  (x  -  xa)/'(x«)  +  (x  f'\q{x))  (85) 

where  q(x)  lies  between  z  and  xa.  If  xa  is  a  close  approximation  to  r,  the  quadratic 
term  can  be  neglected  and  equation  (85),  evaluated  at  x  =  r,  gives 

0  «  /(xa)  +  (r  -  xa)/'(xa)  (86) 


Solving  for  r 


r  «  x„  — 


/(*») 

I'M 


(87) 
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This  sets  the  stage  for  an  iteration  procedure,  in  order  to  determine  the  root,  r,  that 
starts  with  an  initial  approximation  zo  and  generates  the  sequence  x*,  defined  by 

/(x*_i) 

(88) 

Newton’s  method  can  be  generalized  to  a  system  of  m  non-linear  algebraic 
equations  of  the  form  /«(x/)  =  0,  where  #e=l,. .  .,m  and  1=1,. .  .,m,  or,  in  vector  form 
F(X)  =  0,  where  X  is  the  vector  that  contains  the  unknowns.  In  the  present  appli¬ 
cation,  /*  =  0  represent  the  five  equations  of  our  problem,  written  in  finite  difference 
form  at  n  grids  in  the  y  direction  across  the  LLBL,  at  a  particular  x  location,  along 
with  one  equation  at  the  boundary  grid  n-fl;  xj  represent  all  the  variables,  i.e.,  the 
unknown  values  of  u,  u,  B,  p  and  p  at  all  n  grid  points,  and  un+1>  Thus  we  have  m 
=  (5xn)+l  in  the  present  application. 

The  functional  iteration  procedure  evolves  by  first  selecting  X0  and  then  itera¬ 
tively  generating 

** - **-■  -  juhF<x'-'}  (89> 

for  k  >  0.  Here,  J(X)  =  dfi/dxj  is  the  Jacobian  of  F ,  the  detailed  expression  for 
which  is  given  below.  One  can  prove  that  the  sequence  generated  by  this  iteration 
gives  quadratic  convergence  to  the  solution,  P,  of  F(X)  =  0,  provided  that  the 
initial  guess,  Xo,  is  a  sufficiently  close  approximation  to  the  actual  solution  [Burden 
and  Faires ,  1989].  In  order  to  satisfy  this  condition,  we  take  advantage  of  the  slow 
evolution  of  the  LLBL  downstream  by  using  the  values  of  all  variables  in  the  previous 
step  as  an  initial  guess  for  a  new  step  in  the  downstream  direction. 

In  practice,  explicit  computation  of  1  /  J(X)  is  avoided  by  performing  the  oper¬ 
ation  in  a  two-step  manner.  First,  a  vector  Y  =  Xk  -  Xk.x  is  found  which  satisfies 

J(Xk-i)Y  =  -F(Xk-t)  (90) 
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After  this  has  been  accomplished,  the  new  approximation,  Xk,  is  obtained  by  adding 


Y  to 

The  vector  Xk  contains  the  solution  ( [uj ,  m,  Bx ,  pi ,  pi , . . . ,  un,  v„,  Bn ,  ,  pn ,  «n+i  ]t ) r 

after  each  iteration  step,  k.  The  vector  F(X)  contains  the  left-hand  side  of  the 
(5xn)+l  finite  difference  equations  described  earlier.  The  Jacobian  matrix  J(X)  is 
of  the  form 


M  QLl  g/i 

5u7  dvi  5ti„+i 


I  dftn+ 1  dftn+l  dfin+l  I 

L  ~^r  -  -xztr  J 

Since  all  of  the  finite  difference  equations  involve  the  variables  at  no  more  than 
three  neighboring  grid  points  in  the  y  direction,  it  is  possible  to  arrange  the  above 
matrix  in  a  banded  form  with  bandwidth  equal  to  15.  Notice  from  equation  (89)  that, 
at  every  step  of  the  iteration  procedure,  after  the  initial  guess,  X0,  both  J(Xk-i)  and 
F(Xk-i)  are  known,  and  therefore,  a  linear  algebraic  system  has  to  be  solved  for  Xk. 
This  is  done  by  an  IMSL  (Mathematical  and  Statistical  Library)  Library  subroutine, 
that  performs  LU  decomposition  of  the  above  matrix  and  includes  partial  pivoting. 
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Fig.  7.  A  schematic  of  the  computational  domain  and  the  asymptotic  domain 
in  the  equatorial  plane. 
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APPENDIX  4:  BENCH  MARK  TESTS 

The  system  of  equations  (15),  (17),  (18),  (20)  and  (26)  possesses  certain  self- 
similar  solutions.  In  this  Appendix,  these  solutions  are  derived  and  used  as  bench¬ 
marks  for  the  numerical  code  described  in  Appendix  3. 

We  first  examine  the  mass  conservation  law.  The  x  and  y  components  of  the 
velocity,  and  the  plasma  density  are  assumed  to  have  the  following  variation  with  x 
and  y: 

Vro  =  xmfiv),  Vyo  =  xpg(r)),  p0  =  x,r(??),  7  =  (92) 

Here,  /,  g,  and  r  are  functions  of  the  independent  variable  7  only,  and  m,  p ,  q  and 
n  are  exponents  to  be  determined  in  such  a  way  that  all  powers  of  x  are  cancelled  in 
the  system  of  equations  mentioned  above.  A  solution  derived  in  this  manner  has  the 
following  property:  one  may  obtain  the  solution  at  every  x  location  from  the  solution 
at  any  other  x  location  by  (a)  multiplying  by  x  raised  to  the  appropriate  power  and 
(b)  rescaling  the  y  axis  so  that  the  ratio  y/xB  matches  the  similarity  variable,  7.  By 
substituting  (92)  into  the  mass  conservation  law,  equation  (15),  we  get 

(m  +  q)fr  -  n7(/r)'  +  (yr)'  =  0  (93) 

The  differentiation,  denoted  by  the  prime,  is  with  respect  to  the  only  remaining 
independent  variable,  namely,  7.  The  requirement  that  the  powers  of  x  cancel  results 
in  the  following  condition  for  the  exponents: 

m  -  1  =  p  —  n  (94) 

To  derive  an  ordinary  differential  equation  from  the  x-momentum  equation, 
we  assume  that  the  total  pressure  in  the  equatorial  plane,  /^(x),  has  the  following 
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variation: 


/>«,(*)  =  1 2m*’C„  (95) 

where  Cp  is  a  constant,  independent  of  both  x  and  q.  From  the  y-momentum  equation, 
(17),  it  then  follows  that 


Po  =  x2m+qx(q),  Bt 0  =  xm+qt7b(r)) 

and  the  y-momentum  equation  takes  the  form 

b2  _ 

X  +  2^-C'’ 


(96) 


(97) 


Furthermore,  we  assume  that  the  conductivity,  a,  the  current,  jv  and  the  viscosity, 
t)0  have  the  following  variation: 


ff  =  x"E(,),  =  >)),  Vo  =  x‘M(<l) 

With  the  above  assumptions,  the  x- momentum  equation,  (26),  becomes 
mrf 2  -  nqrf  f  +  rgf  =  -(2m  +  q)Cp  -  E62/ 

+  (  S00/00600  +  •/«,+  ^(m  +  |)600^  b  +  M'f'  +  M/" 

The  requirements  on  the  exponents  are: 

2m  +  q  —  1  =  3m  -f  q  +  fi  =  s  +  m  +  q/2  =  a  —  2n  +  m 
from  which  we  get 

q 

fi  =  —  1  —  m,  s  =  m  +  -  —  1,  a  =  m  +  2n-f-g  —  1 

Similarly,  the  isentropic  law,  equation  (18),  becomes 
(2m  +  q  —  7?)/ ~  —  ni)f  (~)  +  9  (£)  =0 


(98) 


(99) 


(100) 


(101) 


(102) 
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(103) 


and  the  induction  law,  equation  (20),  becomes 

(">  - 1)/‘  —  1/  (*)  +»(;)=» 

No  further  restrictions  on  the  exponents  are  imposed  by  the  latter  two  equations. 

The  system  (93),  (97),  (99),  (102),  and  (103)  consists  of  five  ordinary  differential 
equations  for  /,  g,  r,  t  and  b,  with  respect  to  the  only  independent  variable,  q;  it  can 
be  integrated  subject  to  five  boundary  conditions  (one  equation  is  of  second  order 
and  one  is  algebraic).  It  contains  seven  exponents,  namely,  m,  n,  p,  q,  ft,  s  and  a, 
which  are  constrained  by  the  four  conditions  (94)  and  (101),  so  that  three  exponents 
can  be  specified. 

If  the  total  pressure  in  the  magnetosphere  is  a  function  of  x,  the  x  variation  of 
Po,  po  and  BM  is  given  by  (92)  and  (96).  The  solution  that  was  used  for  benchmarking 
of  the  numerical  code  in  the  computational  box  is  of  this  form.  The  exponents  were 
chosen  as  follows:  n  was  chosen  to  be  zero,  which  represents  a  “constant  thickness” 
boundary  layer;  in  this  case  q  =  y  and  no  scaling  in  the  y  direction  is  needed  in  order 
to  relate  the  solutions  at  different  x  locations;  m  was  chosen  equal  to  1,  and  q  was 
chosen  equal  to  2  so  that  the  magnetopause  boundary  is  a  stream  line  with  ^(0)  =  0. 
With  n  as  0,  m  =  1,  and  q  =  2  it  follows  from  (94)  and  (101)  that  p  =  0,  ft  =  —2, 
s=l,  and  a  =  2.  One  can  see  that  with  those  parameters,  equations  (97),  (102),  and 
(103)  require  solutions  for  ir,  r,  and  b  that  are  independent  of  q.  The  five  variables 
then  take  the  form 

v*o  =  xf{y),  Vyo  =  g(y),  p0  =  x2Cr,  po  =  xACr  Bt o  =  x2Cb  (104) 

Notice  from  (98)  that  the  conductivity,  the  current  flowing  into  the  magnetosphere, 
and  the  viscosity  are  functions  of  x.  The  viscosity  function,  M,  was  chosen  to  be  in¬ 
dependent  of  y.  The  resulting  system  of  equations  was  solved  subject  to  the  following 
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boundary  conditions:  /( 0)  =  —180  km/sec,  /( oo)  =  0,  g( 0)  =  0,  £>o(0,y)  =  44  nT, 
and  y)  =  0.62  nPa. 

The  solution  was  obtained  by  integrating  (93)  and  (99);  an  IMSL  (Mathematical 
and  Statistical  Library)  Ordinary  Differential  Equation  solver  routine  was  used  for 
that  purpose.  This  solution  was  then  used  as  the  initial  step  in  the  code  described  in 
Appendix  3,  at  x  =  —10  Re .  Boundary  conditions  and  other  parameters  were  chosen 
according  to  the  requirements  for  the  self-similar  solution  described  above.  Figure 
8  shows  the  profiles,  obtained  by  the  numerical  code,  of  v*o,  *V>,  ^*o>  Po,  and  po, 
versus  y  in  the  computational  box,  at  two  different  x  locations,  namely,  after  3  Re 
and  after  6  Re  from  the  beginning  of  the  calculation.  At  each  location,  the  profile 
from  the  independent  integration  of  the  ODE’s,  scaled  to  the  appropriate  x  location 
is  superimposed.  The  difference  between  the  two  solutions  is  of  the  order  of  one  part 
in  10'*  and  is  indistinguishable  in  the  plots. 

It  should  be  noted  that  this  case  is  not  a  complete  benchmark  of  the  code  because 
po,  po  and  were  chosen  to  be  functions  of  x  only:  those  variables  are  constant 
with  y.  However,  it  can  be  argued  that  this  is  not  a  major  restriction  because  (a) 
the  pressure  balance  in  the  y  direction,  which  is  explicitly  incorporated  in  the  code, 
is  accurately  satisfied  at  every  step  in  the  —  x  direction,  and  (b)  the  validity  of  the 
finite  difference  form  for  the  convective  terms  in  the  induction  law  (the  isentropic  law 
is  expressed  similarly),  when  Bz o  is  non-constant  with  y,  has  been  verified  separately. 

It  is  emphasized  that  in  this  benchmark  test,  an  explicit  boundary  condition 
on  the  value  of  v*o  at  the  magnetospheric  boundary  was  used  in  the  numerical  code, 
instead  of  the  mixed  boundary  condition  derived  by  matching  the  solution  at  the 
computational  box  with  an  asymptotic  solution  that  extends  into  the  magnetosphere. 
This  serves  the  purpose  of  verifying  the  finite  difference  representation  of  the  equations 
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described  in  Appendix  3.  It  was  found  that,  when  a  computational  solution  is  matched 
to  an  asymptotic  solution  at  y  =  y&,  the  overall  accuracy  of  the  code  is  somehow 
influenced  by  the  choice  of  the  boundary  condition  applied  at  y  =  y*.  The  method 
we  chose  was  to  require  a  solution  that  would  guarantee  continuous  first  and  second 
derivatives  of  the  velocity  component  v^q  at  y  =  y&  (in  addition  to  continuity  of  t?*o 
itself),  so  that  the  first  derivative  of  the  current  at  the  top  of  the  boundary  layer 
is  also  continuous,  as  explained  in  Appendix  1.  The  boundary  condition  that  was 
implemented  is  given  by  equation  (33)  and  the  solution  in  the  asymptotic  box  is 
given  by  (29),  where  a(x)  has  been  slightly  adjusted  from  the  value  given  by  (30) 
to  allow  for  a  continuous  second  derivative  of  u*o.  For  the  benchmark  test  of  the 
entire  code,  including  the  asymptotic  box,  a  self-similar  solution  with  constant  total 
pressure  along  the  equatorial  LLBL  was  chosen,  as  described  below. 

If  the  total  pressure  in  the  magnetosphere,  P*,,  is  assumed  to  be  independent 
of  x,  then  the  density,  the  plasma  pressure,  and  the  z  component  of  the  magnetic 
field  should  also  be  independent  of  both  x  and  y  in  order  to  satisfy  equations  (17), 
(18),  and  (20).  The  mass  conservation  law,  which  is  then  of  the  form  =  0,  i.e., 
it  represents  an  incompressible  solution,  becomes 

mf  -  nrjf  +  g'  =  0  (105) 

and  the  x-momentum  equation  becomes 

CT{mp  -  nrjf  f  +  gf)  =  -EC?/  +  (Woo  Cb  +  +  M'f  +  Mf "  (106) 

Here,  Cr  and  C&  denote  the  constant  values  of  p0  and  Bxq.  The  relation  (94)  for  the 
exponents  is  still  valid  and  the  relation  (100)  becomes 

2m-l  =  m  +  /?  =  s  =  a-2n  +  m  (107) 
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The  solution  of  (105)  and  (106)  that  was  used  for  benchmarking  has  an  asymp¬ 
totic  behavior  as  y  — ►  oo  such  that  jy =  0,  u*oo  =  0.  Abo,  the  values  m=l,  n=0 
were  used.  Under  those  assumptions  the  solution  to  (105)  and  (106)  can  be  found 
analytically;  it  is  given  by 

f  =  Ae $  =-(«"*-  1)  (108) 

c 

where  A  is  given  by  the  boundary  condition  on  /  at  y  =  0  and  c  satisfies  the  relation 

A  =  s§f  +  r/  (109) 

Figure  9  shows  the  profiles  of  Vzo  and  Uyo  versus  y,  at  a  distance  of  10  Re  downstream 
from  the  beginning  of  the  calculation.  The  computational  box  ends  at  the  location 
where  v*o  =  -3  km/s.  As  in  Figure  8,  two  profiles  have  been  superimposed,  one 
from  the  numerical  code  and  the  other  from  the  independent  ODE  calculation.  The 
difference  between  the  two  results  is  of  the  order  of  one  part  in  102. 

4.1  Self-Similar  Solutions  of  the  Extended  Model 

The  extended  model  described  in  Appendix  2,  that  includes  field-aligned  po¬ 
tential  drops  in  the  coupling  region,  possesses  self-similar  solutions  as  well. 

In  the  following,  Vzq,  Vyo,  Bx 0,  po,  po,  and  tj0  are  assumed  to  have  the  functional 
behavior  described  earlier  in  this  Appendix.  Furthermore, 

Bxl  =  xed(rj),  jNf  =  xcJ||(»7),  Bi  =  xblBi(rj),  ^  =  xXidx{r]), 

s p  =  x'£(i?),  H0  =  x^HM,  Byl  =  x'byiri),  vtX=xvz{r})  (110) 

With  the  above  assumptions  the  x-momentum  equation,  (16),  becomes 

mrf 2  -  mjr//'  +  rgf  =  -(2m  +  q)Cp  +  - l—bd  +  A /'/'  -(-  M f"  (111) 

tip  o 
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where  the  requirements  on  the  exponents  are 

2m  +  q  —  1  =  m  +  ^-M  =  m--2n  +  a  (112) 

A 

The  y-momentum  equation  remains  as  in  (97),  and  the  isentropic  law  remains  as  in 
(102).  The  requirements  on  the  exponents  from  this  law  are  again 

m  —  1  =  p  —  n  (113) 

The  induction  law,  equation  (19),  becomes 

(2m  +  i)bf  -«,(»/)'  +  (is)'  =  0  (114) 

while  equation  (67)  becomes 

J||  =  s(y*)1  +  (/*)’)  +  y-fa  +  fb)  (115) 

The  requirements  on  the  exponents  are 

a  +  26j  +  2x\  =  2m  4-  2n  +  q,  c  —  2m  +  n  +  q/2  (116) 

The  current  continuity,  equation  (63),  becomes 

along  with  with  the  requirement 

h  —  n  +  e  =  c  —  bi+m  +  q/2  (118) 

The  mass  conservation  law,  equation  (47),  gives 

(q  +  h  +  m)rHnf  -  ( nr)  -  1  ){rH„f)'  =  0  (119) 

and  the  2-momentum  equation,  (50),  becomes 

2M,r°o  -  »)  =  (12°) 
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along  with  the  requirement 


e  +  /i  =  to  +  q/2  (121) 

Equation  (46)  becomes 

yfj 

-  nrjH^)  -  gW,  +  ^  =  0  (122) 

along  with  the  requirement 

m  —  1  =  p  —  n  =  v.  (123) 

Finally,  the  magnetic  flux  conservation  law,  equation  (52),  becomes 
ed  —  nrid '  +  6'  =  0  (124) 

along  with  the  requirement 

t-  n  =  e  —  1.  (125) 


The  algebraic  system  (112),  (113),  (116),  (118),  (121),  (123),  and  (125)  gives 
e  =  m  +  q/2  —  1 
a  =  m  +  2n  +  g  —  1 
p  =  to  +  n  —  1 
c  =  2  to  +  n  +  q/2 
bi  =  2  to  +  2n  +  q/2 
h  =  1 

a  +  2x\  =  —2m  —  2n 
v  =  to  —  1 
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t  =  2m  +  q/ 2  —  2 


(126) 


The  above  nine  relations  constrain  the  thirteen  parameters  that  are  present  in  the 
system  of  ordinary  differential  equations  derived  in  this  section.  Four  of  those  pa¬ 
rameters  are  free,  and  one,  namely  h,  is  required  to  have  a  fixed  value.  Note  that  this 
results  in  a  somewhat  unphysical  variation  with  x  of  the  height  Ho',  nevertheless,  the 
solution  can  be  used  for  the  purpose  of  benchmarking. 

4.2  Other  Benchmarks 

As  already  mentioned  in  Appendix  1,  the  system  has  been  benchmarked  against  the 
one-dimensional  solution  (  which  is  independent  of  x)  derived  by  Sonnerup  [1980], 
under  the  appropriate  assumptions.  First,  the  solution,  which  is  given  by  equation 
(42),  was  used  as  an  upstream  velocity  profile:  it  was  found  to  remain  unchanged, 
for  as  long  a  distance  downstream  as  desired.  Second,  a  different  initial  profile  was 
imposed  and  it  was  found  that,  with  increasing  |x|,  the  system  eventually  relaxed  to 
Sonnerup’s  solution,  with  vv  =  0,  as  required. 
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Fig.  8.  Self-similar  benchmark  solution:  plotted  are  profiles  of  Vxo,  v^o  (top  two 
panels),  pQ,  (middle  two  panels),  and  po  (bottom  panel),  versus  y,  after  a  flow 
distance  of  3  Re  and  6  Re  from  the  upstream  location.  Each  profile  is  calculated  from 
the  numerical  code  and  from  an  ODE  solver  independently,  and  the  two  solutions  are 
superimposed. 


vx  (Jcm/3) 


0  y  1  re 

Fig.  9.  Incompressible  self-similar  solution.  The  two  panels  show  the  profiles 
of  Vxo  and  Vyo,  respectively,  versus  y,  at  flow  distance  of  10  Re  from  the  beginning  of 
the  calculation;  Bz 0,  po,  and  p0  are  constant.  The  solutions  given  by  the  numerical 
code  and  by  independent  integration  of  the  ODE  are  superimposed.  The  computa¬ 
tional  box  ends  at  0.3  Re,  and  an  asymptotic  solution  is  pasted  to  the  computational 
solution  at  that  location. 
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ABSTRACT 

A  mathematical  model  for  a  force-free  boundary  layer  (FFBL)  has  been 
developed  to  calculate  the  deflection  of  the  earth’s  magnetic  field  due  to  ^uasi-steady 
field-aligned  currents.  The  model  may  be  used  to  determine  the  magnetic  field  struc¬ 
ture  and  mapping  between  the  equatorial  magnetosphere  and  the  ionosphere.  The 
geometrical  volume  of  interest  extends  between  two  magnetic  flux  surfaces,  with  a  low 
altitude  boundary  representing  the  ionosphere  and  a  high  altitude  boundary  represent¬ 
ing  the  interface  with  a  model  for  field-aligned  current  generation.  The  mathematical 
formulation  is  general  and  may  be  implemented  numerically  for  any  magnetic 
geometry  for  which  locally  orthogonal  coordinates  can  be  defined.  A  numerical  imple¬ 
mentation  of  the  model  and  its  application  to  dayside  region  1  currents  are  described. 
The  illustrative  application  in  dipole  magnetic  geometry  suggests  that  typically 
observed  dayside  region  1  currents  produce  a  maximum  (upper  limit)  azimuthal 
deflection  of  dipole  field  lines  of  about  26°. 


1  Now  at  Hewlett-Packard  Company,  Workstation  Division,  Chelmsford,  Mas¬ 
sachusetts. 
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1.  Introduction 


The  description  of  magnetic  field-aligned  currents  (FAC)  presents  a  difficult  problem  to 
the  magnetospheric  modeler  because  the  closure  paths  of  the  currents  are  determined  funda¬ 
mentally  by  the  local  plasma  and  fluid  behavior.  In  the  framework  of  one-fluid  magnetohydro¬ 
dynamics  (MHD),  FAC  closure  can  be  effected  by  plasma  polarization  (finite  inertia),  plasma 
diamagnetism  (finite  thermal  energy),  and  plasma  electrical  conductivity  (finite  collision  time). 
The  last  effect  dominates  in  the  ionosphere  while  the  first  two  are  important  in  the  equatorial 
magnetosphere.  In  the  intermediate  altitude  region,  between  the  ionosphere  and  the  equatorial 
magnetosphere,  the  transverse  elements  of  the  plasma  conductivity  tensor,  the  ratio  of  plasma 
pressure  to  magnetic  pressure,  and  the  Alfv6n  Mach  number  are  all  small.  Therefore,  the 
electrical  currents  that  flow  in  the  intermediate  altitude  region  are  force-free  to  a  good  approx¬ 
imation  and  remain  field-aligned  until  reaching  the  ionosphere  or  the  equatorial  magnetos¬ 
pheric  region  where  they  are  diverted  into  perpendicular  currents. 

The  purpose  of  this  paper  is  to  describe  a  new  technique  for  calculating  quasi-steady 
field-aligned  currents  in  the  intermediate  altitude  region,  and  the  magnetic  deflections  pro¬ 
duced  by  them,  when  the  currents  flow  in  a  relatively  thin,  but  finite  thickness,  layer.  Stem 
[1993]  has  recently  described  a  technique  for  calculating  the  perturbing  effects  of  field-aligned 
currents  on  the  magnetic  field  in  regions  outside  the  field-aligned  current  layer,  i.e.,  in  regions 
where  the  field-aligned  current  is  zero.  The  magnetic  field  is  evaluated  by  Stem  assuming  that 
the  FAC  flows  along  the  unperturbed  magnetic  field  lines  in  a  zero  thickness  (5-function) 
sheet  or  shell.  When  the  current  layer  has  a  finite  thickness,  the  FAC  must  flow  along  mag¬ 
netic  field  lines  that  are  determined  by  solving  for  the  FAC  path  and  the  magnetic  field  lines 
simultaneously.  Therefore,  the  model  described  in  this  paper  effectively  resolves  the  5- 
function  sheet  currents  assumed  by  Stem.  Because  the  currents  are  assumed  to  be  entirely 
field-aligned  in  the  thin  layer  of  interest  here,  we  refer  to  the  model  as  a  force-free  boundary 
layer  (FFBL). 

The  utility  of  the  FFBL  model  described  here  is  twofold.  It  can  be  used  to  map  FACs 
outward  from  the  ionosphere  where  statistical  synoptic  data  on  their  spatial  distributions  are 
available.  The  model  may  also  be  used,  in  modular  form,  to  connect  a  high  altitude  magnetos¬ 
pheric  dynamo  region  to  its  ionospheric  load.  For  example,  we  envision  (though  have  not  yet 
accomplished)  connecting  the  FFBL  model  with  the  low-latitude  boundary  layer  model 
developed  recently  by  Drakou  et  al.  [1994].  Other  applications  for  connection  to  nightside  or 
cusp  region  dynamos  are  also  possible. 

In  this  brief  report,  the  basic  assumptions  and  mathematical  formulation  of  the  FFBL 
model  are  described,  along  with  a  particular  numerical  implementation.  Use  of  the  basic  tech¬ 
nique  is  illustrated  by  mapping  the  dayside  region  1  currents  [lijima  and  Potemra,  1976]  out¬ 
ward  from  the  ionosphere  in  an  (oversimplified)  approximation  where  the  magnetic  field  out¬ 
side  the  FFBL  remains  dipolar. 

For  reference,  the  field  outside  the  FFBL  is  identified  in  this  paper  as  the  exterior  field. 
The  exterior  field  occupies  regions  that  are  both  closer  to  and  further  from  earth  than  the 
FFBL.  We  assume  throughout  the  paper  that  a  locally  orthogonal,  ‘magnetic’  coordinate  sys¬ 
tem  exists,  in  which  one  of  the  coordinates  lies  along  the  exterior  magnetic  field,  a  second  is 
normal  to  the  FFBL  surface  (also  a  magnetic  flux  surface),  and  the  third  completes  the  orthog¬ 
onal  set.  In  fact,  one  can  always  identify  a  system  of  locally  orthogonal  unit  vectors  at  any 
point  on  a  magnetic  flux  surface  (e.g.,  the  Frenet-Serret  set  or  the  LMN  set  often  used  in 
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magnetopause  studies)  for  the  puiposes  of  resolving  field  vectors  into  locally  orthogonal  com¬ 
ponents  and  performing  vector  algebra  on  them.  However,  representation  of  vector  differential 
operators  requires  specification  of  the  transformation  equations  that  relate  appropriate  mag¬ 
netic  flux  coordinates  to  Cartesian  coordinates  and  which  determine  the  metric  tensor  for  the 
transformation.  Transformations  that  remain  locally  orthogonal  in  the  magnetic  flux  coordi¬ 
nates  are  known  to  exist  for  plane  fields  and  axisymmetric  fields,  of  which  the  dipole  field  is 
a  familiar  example,  but  locally  orthogonal  transformations  can  be  difficult,  if  not  impossible, 
to  find  for  realistic  magnetic  field  configurations.  As  a  consequence,  we  are  currently  in  the 
process  of  generalizing  the  FFBL  formulation  described  here  and  expect  in  a  future  paper  to 
report  on  a  FFBL  model  based  on  non-axisymmetric,  non-orthogonal  magnetic  flux  coordi¬ 
nates. 


2.  Model 


2.1  Force-Free  Condition 

In  the  region  of  interest,  the  ratio  of  plasma  pressure  to  magnetic  pressure,  is  assumed  to 
be  much  smaller  than  one,  as  is  the  ratio  of  plasma  kinetic  energy  density  to  magnetic  pres¬ 
sure.  Therefore  the  Lorentz  force  term  dominates  the  pressure  gradient  and  inertial  terms  in 
the  MHD  equation  of  motion.  The  momentum  equation  then  reduces  to  the  so-called  force- 
free  condition  [e.g..  Priest,  1982]  in  terms  of  the  plasma  current  density  j  and  the  magnetic 
field  B: 

jxB  =  0.  (1) 

From  the  steady-state  Ampere’s  law, 

Poj  =  VxB  (2) 

where  Po  is  the  permeability  of  free-space.  The  force-free  condition  can  be  alternatively 
expressed  as 

Mol  =  oB  (3) 

where  a(r)  is  a  scalar  function  of  the  spatial  coordinates.  Combining  (2)  and  (3)  gives 

V  x  B  =  aB.  (4) 

The  divergence  of  (4),  together  with  V  •  B  =  0,  shows  that  a(r)  is  constant  along  magnetic 
field  lines: 

B  ■  Va  =  0.  (5) 


2.2  Geometrical  Domain 

After  implementing  the  approximations  described  in  Section  2.3  below,  a  solution  to  (4) 
and  (5)  will  be  developed  in  a  general  curvilinear  coordinate  system  with  the  proviso  that  the 
coordinates  be  locally  orthogonal.  The  curvilinear  coordinates  are  as  follows:  x,  is  the  gen¬ 
eralized  ‘azimuthal’  coordinate;  x2  is  the  flux  surface  label  and  decreases  in  value  when  going 
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from  inner  flus  surfaces  to  outer  flus  surfaces;  x3  varies  along  the  exterior  magnetic  field.  The 
unit  vectors  k2  and  ft3  form  a  mutually  orthogonal  set.  As  mentioned  above,  the  term 
‘exterior’  is  used  to  indicate  the  magnetic  field  outside  and  at  the  edge  of  the  FFBL  where  the 
FAC  is  zero.  In  general,  the  exterior  field  is  also  influenced  by  the  force-free  current;  these 
effects  can  be  modeled  in  the  region  outside  the  force-free  boundary  magnetic  by  calculating 
the  magnetic  field  perturbations  produced  by  a  zero  thickness,  force-free  current  sheet  (or 
shell),  i.e.,  a  5-function  current  sheet.  This  current  sheet  can  then  be  resolved  into  a  finite 
width,  force-free  current  layer  using  the  procedures  described  below.  The  geometrical  region 
of  interest  is  therefore  a  thin,  but  finite  thickness  layer  where  the  force-free  current  is  nonzero. 

This  geometrical  volume  is  illustrated  in  Figure  1.  It  extends  between  two  magnetic  flux 
surfaces  (to  which  the  exterior  field  is  everywhere  tangential),  with  a  low  altitude  boundary 
representing  the  ionosphere  and  a  high  altitude  boundary  representing  the  interface  with  a 
model  for  field-aligned  current  diversion  or  generation  in  the  magnetosphere.  The  inner  and 
outer  magnetic  shells  correspond,  respectively,  to  the  maximum  and  minimum  values  of  the 
coordinate  jc2.  For  applications  of  the  model  to  portions  of  the  dayside  magnetospheric  boun¬ 
dary  layer  region  on  closed  field  lines,  the  outer  magnetic  shell  might  be  taken  as  the  magne¬ 
topause  surface  and  its  projection  along  field  lines  from  the  equatorial  plane  to  the  ionosphere. 
Although  the  numerical  calculations  described  later  in  this  paper  are  based  on  dipole  magnetic 
geometry,  the  mathematical  formulation  is  general  and  may  be  implemented  numerically  for 
any  magnetic  geometry  for  which  locally  orthogonal  coordinates  can  be  defined. 

23  Boundary  Layer  Approximation 


The  method  of  solution  to  (4)  and  (5)  makes  use  of  a  narrow-channel  or  boundary-layer 
approximation  for  which  spatial  changes  across  the  thin  layer  in  the  x2  direction  occur  on  a 
scale  length  5  that  is  much  less  than  the  scale  length  L  for  changes  in  the  x ,  direction  or  the 
scale  length  H  in  the  x3  direction.  The  x ,  and  x3  coordinates  vary  on  surfaces  that  are 
tangential  to  the  thin  layer.  The  boundary  layer  approximation  implies  that 

a/dx,-l/L,  a/&c2~l/5,  d/dx3~l/H  (6) 

where  5  <i  -  H .  Assigning  a  characteristic  value  B0  for  the  field  components  Bx  and  B3 
in  the  ft,  and  ft3  directions,  the  solenoidal  condition  V  B  =  0  implies  the  following  scaling 
relations: 

B,  -  B3  -  B0;  B2-|bo<Bo.  (7) 


B2  is  the  component  of  B  in  the  x2  direction. 

Expressing  the  curl  in  Ampere’s  law  in  a  locally  orthogonal  coordinate  system  with 
metric  scale  factors,  /i(  =  |dr/dxj,  provides  the  following  three  equations: 


1 

H0M3 


72  = 


1 

M  1^3 


h 


1 

m/*  1^*2 


d(h  3B3) 

d(h2B2) 

dx2 

dx3 

d(h  1 B  |) 

d(h3B3) 

dx3 

dxt 

d(h2B2) 

d(h,Bi) 

Bxt 

dx2 

(8) 

(9) 

(10) 


64 


Using  the  boundary  layer  scaling  relations  defined  by  (6)  and  (7),  and  neglecting  terms 
of  O (82/L2)  in  Ampere's  law,  equations  (8)  and  (10)  become: 


.  _  1  ^3g3> 

M2*3  dx2 

■  ~  1  1^1) 

dx2 


(11) 

(12) 


Comparison  of  the  second  equation  (9)  of  Ampere’s  law  with  first  (8)  and  third  (10)  equations 
shows  that 


h 
i  1 


h 

h 


(13) 


The  current  flow  normal  to  the  boundary  layer  is  therefore  small. 


Combining  (11)  and  (12),  and  using  the  force-free  condition,  jiq}  =  otB,  yields  the  fol¬ 
lowing  equation  for  B3: 


_ 1 _ a_ 

hxh2a  dx2 


h  j  d(h  38  3) 

h2h3a  dx2 


+83  —  0. 


(14) 


An  analogous  equation  can  be  derived  for  B ,.  We  now  further  assume  the  layer  is 
sufficiently  thin  so  that  the  metric  scale  factors  hx  and  h3  are  very  nearly  constant  on  the 
scale  size  for  variations  in  B  across  the  layer:  |82ln/i  |  3| «  |32lnB  I  3|.  With  this  approxima¬ 
tion  (14)  becomes 
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(15) 


2.4  Boundary  Layer  Solution 


To  solve  (15)  a  new  variable,  dx  =  ah2jdx2,  is  defined.  Equation  (15)  then  becomes 

at2B3(x)+B3(T)  =  0, 

and  B  3  is  determined  from  its  solution  as 


*i(t)  = 


dB3(x) 

ax 


(16) 

(17) 


The  general  solution,  with  explicit  dependence  on  xx,  x2,  and  x3  retained,  is 


B  \  =Bw(xi,x3)cos[v(x1,x3)+t(j:,,X2,x3)]  (18) 

B3  =  B00(x1,x3)sin[y(x1,x3)+T(x1,x2,x3)]  (19) 


where 

fl«(*i.*3)  =  yjBl2(x\,x2'0,x3)+B32{xlx2„,x3)  (20) 

y(x,,x3)  =  arctanlB^Xj.x^x^/B^x^x^Xj)]  (21) 

*2- 

x(x,,x2,x3)=  J  a (xx,x2,x3)h2(xx,x2,x3)dx2  (22) 

*2 
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The  constants  of  integration,  B„  and  \y,  are  determined  by  boundary  conditions  on  B  (  and  B3 
on  a  reference  flux  surface  identified  as  x2  =  x2oo.  The  x2  dependence  in  the  metric  scale  fac¬ 
tor  h2  appearing  inside  the  integral  in  (22)  has  been  formally  retained,  although  no  accuracy 
is  actually  gained  by  including  this  dependence.  The  approximation  used  in  going  from  equa¬ 
tion  (14)  to  (15),  in  fact,  allows  one  to  treat  h2  as  constant  across  the  layer  so  that  it  is 
equally  correct  to  set  h2(xx,x2,x3)  ~  h2(xx,x2oo,x3)  in  (22)  and  to  take  h2  outside  the 
integral. 


B2  is  obtained  by  substituting  B{  and  fl3  from  (18)-(22)  in  V •  B  =  0.  After  integration 
with  respect  to  x2 ,  and  using  B2(xl,x2oo,x2)  =  0,  which  must  be  true  if  x2  =  x2mt  is  a  flux  sur¬ 
face,  we  find 


b2 
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(23) 


Note  that  the  metric  scale  factors  may  vary  as  rapidly  (or  slowly)  as  B  in  the  k,  and  ft3 
directions,  so  they  may  not  be  passed  through  the  derivatives  in  the  integrand  in  (23).  It  is 
also  noted  that  hl23  are  in  general  functions  of  xx,  x2,  and  jc3.  However,  as  discussed  above, 
we  gain  no  accuracy  in  retaining  the  x2  dependence  in  the  scale  factors;  they  may  be 
evaluated  at  x{  =  x2oc  and  treated  as  constants  in  performing  the  integration  over  x2.  The 
distribution  of  a(x,,x2,x3)  throughout  the  spatial  region  of  interest  is  not  known  a  priori 
except  on  some  bounding  surface,  say  x3  =  x  j0)  where  xj0)  is  a  constant  value  of  the  coordi¬ 
nate  x  3.  The  volume  distribution  of  a  can  be  determined  self-consistently  along  with  the 
magnetic  field  by  simultaneously  solving  (5)  and  (l8)-(23). 


2.5  Boundary  Conditions 


The  formulation  of  the  force-free  boundary  layer  model  given  above  is  general  and  may 
be  used  when  (i)  the  region  of  interest  can  be  characterized  by  the  force-free  condition  (1)  and 
(ii)  the  field-aligned  current  into  the  region  flows  in  a  thin  layer  so  that  the  boundary  layer 
approximations  (6)  and  (7)  are  appropriate.  Given  these  constraints,  a  unique  solution  to  the 
set  of  equations  (5)  and  (18)-(23)  requires  two  boundary  conditions: 

(A)  specification  of  the  magnetic  field  B(x1,x2w,,x3)  on  a  reference  magnetic  flux  surface 
x2  =  *2-,  and 

(B)  specification  of  a(x1,x2,xj0))  =  \IqJIB  on  a  bounding  surface,  x3  =  x|0),  normal  to  the 
magnetic  field  (a  ‘magnetic  normal  surface’). 

Technically,  a  unique  solution  also  requires  specifying  the  partial  derivative  d(hxh2B3)ldx3  in 
(23)  on  the  bounding  surface  x3  =  x^0) .  We  have  found  in  practice,  however,  that  the  value  of 
this  derivative  on  the  bounding  surface  has  little  influence  on  the  accuracy  of  the  boundary 
layer  solution  (cf.  discussion  in  Sec.  2.6). 

For  applications  to  the  magnetosphere,  we  will  take  the  reference  magnetic  flux  surface 
to  be  located  at  the  inner  edge  of  the  force-free  boundary  layer.  In  dipole  coordinates,  this 
surface  would  correspond  to  the  innermost  L-shell  of  the  FFBL.  While  the  magnetic  field 
earthward  of  the  reference  flux  surface  is  presumed  to  be  known,  this  internal  field,  in  general, 
depends  on  the  currents  that  flow  in  the  FFBL.  Because  the  FFBL  is  a  thin  layer,  its 
influence  on  the  internal  field  could  be  calculated,  approximately,  in  the  limit  where  the  FFBL 
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is  represented  as  a  zero  thickness  (5- function)  shell.  An  internal  held  model,  including  the 
perturbing  magnetic  held  due  held-aligned  currents,  approximated  as  5-function  current  sheets, 
has  been  discussed  by  Stem  [1993].  The  internal  held,  calculated  in  this  way,  then  provides 
boundary  condition  (A)  stated  above;  the  5-function,  held-aligned  current  shell  can  be 
resolved  using  the  FFBL  formulation.  If  better  accuracy  is  desired,  an  iterative  process  must 
be  used.  In  the  illustrative  application  described  in  Sec.  4,  we  do  not  attempt  to  construct  a 
fully  self-consistent  magnetic  held  model.  Instead  the  basic  procedure  for  obtaining  FFBL 
solutions  is  illustrated  using  the  familiar  dipole  held  model  to  specify  boundary  condition  (A) 
and  to  represent  the  locally  orthogonal  curvilinear  coordinate  system  required  in  the  formula¬ 
tion. 

Specification  of  boundary  condition  (B)  requires  either  a  model  for  the  generation  of  the 
held-aligned  currents  or  synoptic  observations  of  the  perturbed  and  background  magnetic  held 
on  a  magnetic  'normal’  surface,  on  which  j,  B,  and  therefore  a  can  be  inferred.  Synoptic 
observations  of  the  currents  are  available  only  at  low  altitudes,  essentially  within  the  iono¬ 
sphere.  For  the  latter  case,  the  FFBL  model,  coupled  with  an  appropriate  internal  magnetic 
held  model,  can  be  used  to  follow  the  observed  low  altitude  currents  to  their  points  of  origin 
in  the  outer  magnetosphere.  Alternatively,  coupling  the  FFBL  model  to  a  model  for  current 
generation,  for  example,  the  low-latitude  boundary  layer  and  region  1  current  generation 
model  developed  by  Drakou  et  al.  [1994],  would  provide  a  means  of  locating  the  ionospheric 
signatures  of  the  magnetospheric  current  generator. 

2.6  Numerical  Algorithm 

The  numerical  algorithm  used  in  our  calculations  assumes  that  the  mapping  between 
locally  orthogonal  curvilinear  coordinates  (xI,x2,x3)  and  a  geometrical  coordinate  system,  for 
example,  Cartesian  (x,y,z)  or  spherical  polar  (r, 8,<j>)  coordinates,  is  known.  In  addition,  the 
metric  scale  factors  hx,h2,h 3  for  the  curvilinear  system  are  assumed  to  be  known.  The  boun¬ 
dary  conditions  stated  above  must  also  be  supplied.  Although  the  surfaces  on  which  the 
boundary  conditions  apply  may  be  quite  complicated  in  the  physical/geometrical  domain, 
these  surfaces  are  simple  planes  in  the  rectangular  computational  domain  spanned  by  the  three 
curvilinear  coordinates.  The  algorithm  solves  for  a{xl,x2,x3)  and  B(xi,x2,x3)  throughout  the 
three-dimensional  FFBL  on  successive  computational  planes,  x3  =  x^k\  where  k  =  0,1,2,.... 
Special  measures  must  be  taken  to  start  the  algorithm  on  the  first  plane. 

On  the  first  plane,  where  a(x1,x2,x|0))  is  known  from  boundary  condition  (B),  we  use 
(18)-(22)  to  calculate  Bx  and  B3  at  grid  points  (x^  ,x^)  where  (ij)  are  integers 
corresponding  to  grid  locations  for  coordinates  xx  and  x2,  respectively.  To  calculate  B2  we 
need  to  evaluate  the  partial  derivatives  and  integral  in  (23).  Initially,  however,  we  have  no 
information  on  the  variation  of  B  in  the  jc3  direction,  which  is  required  to  perform  the  deriva¬ 
tive  d(hxh2B3)/dx3.  To  start  the  algorithm,  we  simply  set  fi2  =  0  on  the  starting  plane. 
Because  B 13  >■  B2  in  the  boundary  layer  approximation,  the  particular  choice  for  B2  on  the 
first  plane  does  not  significantly  influence  the  magnetic  field  mapping;  we  have  found  from 
experience  with  other  starting  procedures  that  the  accuracy  of  the  resulting  magnetic  field 
mapping  differs  from  the  one  used  here  by  terms  of  the  order  of  those  neglected  in  the  boun¬ 
dary  layer  approximation. 
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Our  particular  numerical  method  for  the  rest  of  the  domain  proceeds  as  follows: 

•  Using  (22),  the  metric  scale  factor  h2,  and  the  trapezoidal  rule  of  integration, 
x(xltx2,x^k))  is  calculated  at  every  grid  point  on  the  current  x3  plane,  i.e.  the  x3  =  xjA) 
plane.  Using  (18)-(22),  B ,  and  fi3,  are  then  calculated  at  these  grid  points. 

•  To  calculate  the  B2  component,  we  first  calculate  d(ft2/i3fl,)/dx,  using  central 
differencing  and  d(/>t/t2£3)/dx3  using  forward  differencing.  Then  using  the  trapezoidal 
rule,  we  integrate  (23)  to  obtain  B2  at  grid  points  on  the  x\k)  plane. 

•  The  next  plane,  the  x3  =  xj*  +  l)  plane,  is  chosen  so  that  the  linear  projection  along  B 
from  any  node  in  the  x  plane  to  a  point  in  the  x|* +  plane  stays  within  a  small  pre¬ 
specified  distance  e  of  the  node  in  the  current  x  plane. 

•  Using  the  property  of  force-free  field  lines,  i.e.  a  is  constant  along  the  field  lines,  we  can 
now  determine  values  of  a  at  non-grid  locations  in  the  x3  =  xj*  + 1)  plane  by  following 
the  linearly  projected  field  lines  from  nodes  in  the  x3  =  x3(k)  plane.  We  now  interpolate 
these  now  nonuniformly  spaced  values  of  a  onto  (xj^.x^)  grid  points  in  the  x^*  +  l) 
plane.  To  do  this,  a  is  first  interpolated  onto  an  x2  line  using  a  simple  2  point  interpola¬ 
tion.  Then  a  cubic  spline  interpolation  scheme  is  used  to  interpolate  a  onto  the  grid 
points. 

•  Having  calculated  a  on  all  grid  points,  the  procedure  is  repeated  until  the  final  x3  surface 
of  interest  is  reached. 

2.7  Numerical  Accuracy 

The  accuracy  of  the  numerical  model  was  tested  using  a  benchmark  based  on  a  simple 
analytical  model  of  polar  line  currents  feeding  an  axisymmetric,  5-function  current  sheet  lying 
on  a  magnetic  dipole  L-shell  [Vatan,  1993].  Because  the  polar  line  currents  are  chosen  to 
flow  into  the  south  pole  and  out  of  the  north  pole,  it  can  be  shown  that  the  internal  dipole 
field  remains  exactly  dipolar  so  that  an  exact  analytic  solution  can  be  derived  for  the  purpose 
of  verifying  the  numerical  algorithm.  The  quantity  that  was  compared  is  the  deflection  angle, 
i.e.,  the  azimuthal  deviation  from  ‘orange  segment'  mapping  between  two  surfaces  that  are 
locally  normal  to  the  magnetic  field.  The  accuracy  of  the  above  numerical  procedures  was 
found  to  depend  on  the  following  factors: 

•  Boundary  layer  thickness.  Since  the  model  is  based  on  a  boundary  layer  approximation, 
the  metric  scale  factors  should  not  vary  significantly  across  the  layer.  From  the  bench¬ 
mark  results,  it  was  determined  that  the  variation  in  dipole  scale  factors  across  a  2  RE 
layer  (referenced  to  the  equatorial  plane),  and  straddling  L  =  10,  could  be  as  large  as 
30  %  with  a  corresponding  variation  in  the  deflection  angle  of  9  %  or  less. 

•  Linear  projection  of  magnetic  vectors.  The  maximum  (linearly)  projected  distance 
between  adjacent  magnetic  normal  grid  surfaces  is  controlled  b  preselected  tolerance 
factor  e.  Arbitrarily  small  errors  can  be  achieved  by  dec  a  the  tolerance  factor 
which  decreases  this  distance.  Continuous  field  line  tracing  would  be  realized  in  the 
limit  where  the  tolerance  factor  approaches  zero. 

•  Discretization.  The  error  in  the  deflection  angle  is  very  sensitive  to  changes  in  the  grid 
size  in  the  x3  direction  (distance  between  magnetic  normal  surfaces),  which,  in  turn,  is 
controlled  by  the  above  mentioned  tolerance  factor  e.  The  grid  size  in  the  xt  and  x2 
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directions  ideally  should  be  much  smaller  that  the  gradient  scale  lengths  of  the  initial 
a(r)  distribution. 

•  Differencing  and  integration  schemes.  The  particular  choice  of  differencing  and  integra¬ 
tion  schemes  for  solving  the  force-free  boundary  layer  equations  seems  to  have  a 
minimal  effect  on  the  accuracy  of  the  calculated  deflection  angle.  For  example,  for  the  a 
distributions  considered  in  the  benchmark  study,  there  was  very  little  change  in  the 
results  when  the  accuracy  of  the  integration  scheme  was  improved. 


3.  Application  to  Dayside  Region  1  Currents 

Application  of  the  numerical  algorithm  described  in  the  previous  section  is  illustrated  by 
considering  a  simple  model  for  dayside  field-aligned  currents  observed  at  low  altitudes.  These 
currents  will  be  mapped  outward  to  the  magnetospheric  equatorial  plane  under  the  assumption 
that  the  currents  remain  approximately  force-free  over  a  substantial  length  along  the  magnetic 
field.  For  the  purpose  of  illustration  and  simplicity,  a  dipole  model  is  used  to  specify  boun¬ 
dary  condition  (A)  in  Sec.  2.5,  and  to  represent  a  locally  orthogonal  coordinate  system  based 
on  the  magnetic  geometry.  Hie  familiar  dipole  coordinates  are  specified  in  terms  of  the 
spherical  polar  coordinates  (r,0,0):  *,=<(>,  the  ordinary  azimuth  angle  in  the  x-y  plane; 
x2  =  v  =  sin 20/r,  which  is  constant  along  a  dipole  field  line  and  increases  when  going  to 
lower  L-shells  (note  that  the  L-sheU  is  defined  as  L  -req!RE  where  req  is  the  radial  distance 
to  the  point  where  the  field  line  crosses  the  dipole  equatorial  plane);  and  x3  =  p.  =  cos0/r2, 
which  varies  along  a  dipole  field  line  and  is  constant  along  a  magnetic  potential  orthogonal  to 
a  dipole  field  line.  The  metric  scale  factors  for  dipole  geometry  are 

Av  =  ~jj^-(l  +3cos20)",/i  =  rsind  h ^  =  hj =  — 

A/  is  the  earth’s  magnetic  dipole  moment,  diiu  B  is  the  field  strength. 

According  to  Iijima  and  Potemra  [1976],  large-scale  field-aligned  currents  are  a  statisti¬ 
cally  permanent  feature  of  the  high  latitude  ionosphere.  These  currents  are  concentrated  in  two 
adjacent  annuli  surrounding  the  geomagnetic  pole.  The  higher  latitude  annulus  is  referred  to  as 
the  region  1  current  system;  the  lower  latitude  annulus  is  called  the  region  2  current  system. 
The  region  1  currents  flow  into  the  ionosphere  in  the  morning  sector  and  away  from  the  iono¬ 
sphere  in  the  evening  sector;  the  region  2  currents  flow  in  the  opposite  direction  at  any  given 
local  time.  The  areas  near  local  noon  and  midnight  are  not  so  simply  classified  and  may  be 
strongly  time-dependent.  We  consider  here  the  dayside  region  1  current  system  at  magnetic 
local  times  (MLT)  between  dawn  and  dusk,  excluding  the  region  from  1100-1300  MLT.  Fig. 
2  shows  the  MLT  distribution  of  currents  reported  by  Iijimi  and  Potemra  [1978].  The  upper 
panel  shows  the  statistical  distribution  during  relatively  active  periods  of  geomagnetic  activity; 
the  lower  panel  is  for  quieter  periods. 

The  statistical  local  time  dependence  (i.e.  the  0  dependence)  in  the  postnoon  sector  in 
Fig.  2  ahs  been  fitted  to  a  simple  polynomial  function.  The  polynomial  fit  for  both  strongly 
and  weakly  disturbed  periods,  as  well  as  selected  data  points  from  Fig.  2,  are  shown  in  Fig. 
3.  The  corresponding  curves  for  the  prenoon  sector  shown  in  Fig.  3  are  mirror  reflections  of 
those  in  the  postnoon  sector.  (We  do  not  mean  to  imply  that  the  data  in  Fig.  2  are  also  mirror 
symmetric  about  noon!)  The  figure  includes  fitted  points  between  1100-1300  MLT,  but  these 
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currents  near  noon  are  not  actually  considered  in  the  calculations  described  below.  The 
region  1  field-aligned  currents  flow  in  opposite  directions  on  either  side  of  noon  (as  indicated 
in  Fig.  2)  whereas  Fig.  3  plots  only  the  magnitude  of  the  current. 

For  the  variation  in  field-aligned  current  in  the  V  direction,  we  model  typical  current 
profiles  (actually  magnetic  deflections),  examples  of  which  can  be  found  in  the  paper  by 
Iijima  and  Potemra  [1976]  and  in  many  other  papers  describing  low  altitude  magnetic  field 
data.  Because  satellite  trajectories  are  rarely,  if  ever,  exactly  parallel  to  the  9  direction,  we  do 
not  actually  have  precise  information  on  the  v  variation  of  the  currents  at  low  altitudes.  As  a 
consequence,  we  use  the  data  only  as  a  guide  to  generate  a  current  profile  in  v  that  captures 
the  basic  features  of  the  variation.  The  current  variation  in  v  was  modeled  by  one-half  period 
of  a  sine  function  that  is  zero  at  both  edges  of  the  region  1  current  annulus  and  maximum 
midway  though  it.  The  width  in  v  of  the  annulus  is  chosen  to  correspond  to  1°  in  dipole  mag¬ 
netic  latitude. 

The  product  of  the  current  profiles  in  v  and  0  defines  our  model  for  the  field-aligned 
current  at  the  ionosphere.  The  function  a(0,  v)  on  the  ionospheric  ‘surface’  (taken  to  be,  with 
a  small  error,  a  constant  p  surface  near  1 RE  geocentric)  is  constructed  by  multiplying  the 
field-aligned  current  model  by  Pq  (permeability  of  free  space)  and  dividing  the  product  by  the 
dipole  magnetic  intensity  at  the  earth’s  surface  at  the  appropriate  magnetic  latitude.  This  a 
function  is  given  to  the  algorithm  described  in  Sec.  2  as  input  at  the  ionosphere,  and  the  map¬ 
ping  of  selected  points  on  the  ionospheric  ‘surface’  to  the  magnetospheric  equatorial  plane  due 
to  this  current  is  calculated.  Note  that  we  are  not  physically  connecting  the  FFBL  to  its  mag¬ 
netospheric  dynamo.  The  magnetic  deflection  inferred  from  the  force-free  mapping  of  the 
currents  all  the  way  to  the  magnetospheric  equator  should  therefore  be  regarded  only  as  an 
upper  limit. 

Magnetic  field  maps  at  the  equatorial  plane  for  weakly  disturbed  periods  and  strongly 
disturbed  periods  are  shown  in  Fig.  4  and  5,  respectively.  The  layer  thickness  is  chosen  to  be 
1°  in  dipole  magnetic  latitude,  mid  the  inner  shell  has  been  placed  at  the  magnetic  L-shell, 
L  -  10,  which  locates  the  outer  shell  at  L= 1 1.14.  These  maps  are  displayed  in  a  format  simi¬ 
lar  to  that  use  by  Fairfield  [1991]  and  are  constructed  by  identifying  at  the  ionosphere  (actu¬ 
ally  at  the  earth’s  surface)  ten  different  meridians  of  MLT  from  0700-1100  and  1300-1700  in 
one  hour  steps.  (Only  five  of  the  mapped  meridians  in  Fig.  4  and  3  are  distinct  due  to  the 
assumed  symmetry  about  1200  MLT.)  Force-free  field  lines  intersecting  each  selected  MLT 
meridian  at  1  RE  are  followed  outward  until  they  intersect  the  dipole  equatorial  plane.  In  this 
way,  the  set  of  points  defining  a  MLT  meridian  at  the  earth’s  surface  is  mapped  along  field 
lines  to  the  equatorial  plane.  The  figures  show  the  ten  mapped  MLT  meridians  between 
L  =  10  and  L  =  11.14.  Had  the  field  remained  purely  dipolar  in  the  FFBL,  the  mapped  meri¬ 
dians  would  appear  as  segments  of  radial  spokes,  rather  than  curved  lines,  emanating  from  the 
origin  at  the  indicated  MLTs.  Note  that  a  nonlinear  scale  has  been  used  to  magnify  the  FFBL 
region  between  10  £  L  £  11.14. 

In  the  morning  sector,  the  force-free  current  is  positive  and  into  the  ionosphere,  and 
therefore,  the  deflection  of  the  field  is  toward  dawn;  in  the  evening  sector,  the  current  is  nega¬ 
tive  and  away  from  the  ionosphere,  and  the  deflection  is  toward  dusk.  For  a  specific  example, 
consider  in  Fig.  4  the  curve  labeled  0700  in  the  equatorial  plane,  which  has  been  mapped 
through  the  layer  from  1  RE  altitude.  The  mapped  meridian  in  the  equatorial  plane  lies  further 
away  from  local  noon  than  an  0700  MLT  radial  spoke  would  lie,  at  all  mapped  points,  except 
on  the  inner  L-shell  where  boundary  condition  (A)  of  Sec.  2.5  requires  exact  dipole  mapping. 
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4.  Conclusions 


Given  a  distribution  of  the  ratio  of  the  field-aligned  current  to  magnetic  field  on  a  surface 
whose  normal  is  locally  parallel  B  (a  magnetic  normal  surface),  the  force-free  boundary  layer 
model  and  algorithm  described  here  may  be  used  to  calculate  magnetic  field  deformations  due 
to  force-free  currents  that  flow  in  thin  boundary  layer  regions.  The  model  may  also  be  used 
to  determine  the  magnetic  field  mapping  between  magnetic  normal  surfaces  in  such  regions. 
This  model  is  general  and  can  be  applied  in  any  locally  orthogonal  coordinate  system,  defined 
in  terms  of  the  exterior  magnetic  field  (exterior  means  in  the  region  where  the  force-free 
currents  are  zero,  i.e.,  outside  the  FFBL). 

An  illustrative  application  of  the  model  in  dipole  geometry  to  the  statistical  distribution 
of  region  1  field-aligned  currents  suggests  that  the  average  dayside  region  1  currents  observed 
during  geomagneticaily  active  periods  produce  a  maximum  azimuthal  deflection  of  the  mag¬ 
netic  field  of  26°  when  a  magnetic  field  line  is  followed  from  the  ionosphere  to  its  equatorial 
crossing  point.  The  corresponding  deflection  during  weakly  disturbed  periods  is  about  21°. 
The  maximum  deflections  occur  on  field  lines  where  the  field-aligned  current  density  observed 
at  low  altitudes  maximizes.  These  estimates  of  the  deflection  angle  should  be  regarded  only 
as  upper  limits,  however,  because  the  field-aligned  current  is  not  actually  force-free  along  the 
entire  length  of  the  field  line  from  the  ionosphere  to  equatorial  plane.  To  put  these  estimates 
in  perspective,  the  T87  geomagnetic  field  model  [Tsyganenko,  1987]  exhibits  a  deflection 
angle  of  about  10°  for  a  field  line  that  passes  through  1000  MLT  and  75°  invariant  latitude  in 
the  ionosphere  [Fairfield,  1991].  Although  the  T87  model  does  include  some  effects  of  field- 
aligned  currents,  these  effects  are  introduced  through  an  empirical  polynomial  fitting  function 
in  the  model,  so  it  is  not  known  how  much  of  the  deflection  is  actually  due  to  dayside  region 
1  currents  or  other  magnetic  deformations  cursed,  for  example,  by  the  Chapman-Ferraro 
currents. 

In  considering  future  applications  of  this  model,  two  improvements  are  suggested: 

•  Finding  coordinate  transformations  that  generate  locally  orthogonal,  magnetic  flux  coor¬ 
dinates  based  on  the  exterior  magnetic  geometry  may  not  be  practical  for  realistic  mag- 
netospheric  magnetic  fields;  the  transformations  may  not  even  exist  in  general,  especially 
for  nonaxisymmetric  magnetic  fields.  Formulation  of  the  FFBL  model  in  terms  of 
nonorthogonal  magnetic  flux  coordinates  would  allow  grafting  a  FFBL  model  for  field- 
aligned  currents  onto  more  realistic  magnetic  field  models  such  as  the  semiempirical  T87 
model,  semianalytical  models  [Voigt,  1981;  Hilmer  and  Voigt,  1993;  Schulz  and  McNab, 
1987],  or  gridded  numerical  models  [Toffoletto  et  al.,  1994].  Applications  using  a  more 
realistic  field  model  would  allow  examination  of  the  effects  of  a  magnetic  minimum 
region  near  the  magnetic  cusps  where  the  magnetic  deformations  due  to  field-aligned 
currents  are  likely  to  be  large.  We  are  currently  in  the  process  of  implementing  the  FFBL 
model  in  a  representation  based  on  nonaxisymmetric,  nonorthogonal  magnetic  flux  coor¬ 
dinates. 

•  Although  the  force-free  boundary  layer  model  provides  a  well-defined  procedure  for 
mapping  field-aligned  currents  and  field  lines  between  magnetic  normal  surfaces,  it  does 
not  properly  model  the  dynamo  region  in  the  outer  magnetosphere.  A  better  understand¬ 
ing  of  the  field  line  mapping,  for  example,  between  the  low-latitude  boundary  layer  and 
the  ionosphere,  the  magnetotail  and  the  ionosphere,  or  the  cusp  region  and  the  iono¬ 
sphere  will  require  coupling  the  force-free  boundary  layer  model  to  appropriate  physical 
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dynamo  models.  A  boundary  layer  model  for  the  generation  of  the  region  1  currents, 
based  on  one-fluid  magnetohydrodynamics,  has  been  developed  by  Drakou  et  al.  [1994], 
and  it  is  expected  that  the  two  models  will  be  coupled  in  the  near  future. 
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Figure  2.  Statistical  pattern  of  observed  current  densities  during  (a)  active  periods  and 
(b)  weakly  disturbed  periods  from  Iijima  and  Potemra  [1978]. 
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Data 


Polynomial  fit 


Figure  3.  Polynomial  fit  to  average  'Region  1’  current  density  versus  magnetic  local 
time  for  (a)  active  periods  and  (b)  weakly  disturbed  periods. 
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Figure  4.  Magnetic  field  map  for  the  a  distribution  chosen  to  model  active  periods. 


Figure  5.  Magnetic  field  map  for  the  a  distribution  chosen  to  model  weakly  disturbed 
periods. 
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APPENDIX  6: 

Dynamics  of  Shear  Velocity  Layer  with  Bent  Magnetic  Field 

Lines 

V.  L.  Galinsky  and  B.  U.  0.  Sonnerup 
Thayer  School  of  Engineering,  Dartmouth  College 
8000  Cummings  Hall,  Hanover,  NH  03755-8000 


Abstract.  A  fully  three-dimensional,  magnetohydrodynamic  simulation  of 
velocity-sheared  plasma  flow  in  an  ambient  transverse  magnetic  field  with  bent 
magnetic  field  lines  has  been  performed.  “Ionospheric-like”  boundary  conditions  have 
been  used  for  closing  field-aligned  currents,  the  two  ionospheres  being  represented  by 
conducting  plates  with  constant  resistivity.  We  have  found  a  significant  difference  in  the 
development  of  the  Kelvin- Helmholtz  instability,  compared  to  the  standard  plane  2D 
case  with  a  uniform  transverse  magnetic  field:  the  growth  rate  of  the  instability  drops 
significantly  as  bending  increases.  It  seems  likely  that,  under  conditions  representative 
of  the  Earth’s  low  latitude  boundary  layer  (LLBL),  the  Kelvin-Helmholtz  instability  can 
be  suppressed  completely  by  the  magnetic  field-lin?  tension  if  bending  of  the  magnetic 
field  lines  is  sufficiently  strong.  We  have  also  found  that  a  combination  of  the  tearing 
mode  instability  and  the  Kelvin-Helmholtz  instability  leads  to  the  formation  of  localized 
3D  current/ vortex  tubes,  the  ionospheric  foot  prints  of  which  can  be  considered  as 
possible  models  of  the  auroral  bright  spots  observed  by  the  Viking  satellite.  Quantitative 
comparison  of  our  results  with  satellite  observations  of  velocity-sheared  plasma  flow 
in  the  LLBL  and  its  ionospheric  foot  print  indicates  good  agreement  with  the  chosen 
model  parameters. 

Submitted  to  Geophys.  Res.  Lett.,  April,  1994 
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Introduction 


The  interaction  of  the  supersonic  solar-wind  flow  with  the  Earth’s  magnetosphere 
creates  interfaces  or  narrow  layers  where  either  a  gradual  or  an  abrupt  transition  in 
plasma  and  magnetic  field  properties  occurs  from  interplanetary  values  to  those  of  the 
magnetosphere.  Several  theoretical  steady-state  models  of  one  of  these  regions,  namely 
the  magnetospheric  low  latitude  boundary  layer  (LLBL)  on  closed  field  lines,  have 
been  discussed  in  the  literature  [ Sonnerup ,  1980;  Lotko  et  al.,  1987;  Phan  et  al.,  1989; 
Drakou  et  al .,  1994].  The  steady-state  nature  of  these  models  may  impose  restrictions 
on  their  applicability:  in  particular,  shear  flows  are  subject  to  the  Kelvin-Helmholtz 
(KH)  instability  so  that  the  question  of  the  stability  of  the  LLBL  needs  to  be  addressed. 
Most  results  of  linear  stability  analyses  of  velocity  shear  layers  are  valid  only  when  the 
unperturbed  magnetic  field  lines  are  straight,  and  therefore,  are  not  directly  applicable 
to  the  above  models  where  the  geometry  is  more  complex,  having  parabolic  field  lines 
with  different  curvature  in  different  parts  of  the  layer.  Qualitatively,  such  geometries 
should  have  stability  properties  intermediate  between  the  two  cases  of  magnetic 
field  perpendicular  and  parallel  to  the  plane  of  the  unperturbed  shear  layer.  In  the 
perpendicular  case,  instability  occurs  for  k6  <  2,  k  being  the  wave  number  and  6  the 
shear  width;  in  the  parallel  case,  the  wave  modes  are  stable  for  velocity  jumps  less  than 
twice  the  Alfven  speed  [e.g.,  Miura  and  Pritchett,  1982]. 

Processes  taking  place  in  the  LLBL  are  also  very  important  for  understanding 
of  the  physics  of  mass,  momentum  and  energy  transfer  from  the  solar  wind  to  the 
magnetosphere.  Three  possible  mechanisms  have  been  proposed,  namely  “viscous 
like”  interaction,  magnetic  field  reconnection  and  impulsive  plasma  penetration.  The 
discovery  of  flux  transfer  events  (FTEs)  increased  the  interest  in  the  second  mechanism 
and  stimulated  development  of  new  reconnection  models.  Later  a  combination  of  viscous 
interaction  and  reconnection  has  been  proposed  as  an  expiation  of  FTE  formation  [La 
Belle-Hamer  et  al.,  1988;  Belmont  and  Chanteur,  1989]. 


In  this  paper,  we  investigate  the  influence  of  magnetic  field-line  bending  on  the 
stability  of  a  velocity  shear  layer  by  use  of  full  three-dimensional  MHD  simulations  with 
ionospheric-like  current  closure  boundary  conditions.  Several  examples  of  unstable  as 
well  as  nearly  stable  configurations  will  be  shown  and  the  parameter  range  where  the 
KH  mode  is  strongly  suppressed  by  the  magnetic  field-line  bending  will  be  identified. 
The  formation  of  localized  current/vortex  tubes,  which  provides  a  possible  explanation 
of  FTE-like  signatures,  will  also  be  shown.  Although  there  are  several  restrictions  on 
the  applicability  of  this  model,  we  believe  it  can  be  used  as  a  first  step  in  understanding 
the  properties  of  the  KH  instability  in  a  magnetic-field  geometry  that  is  more  realistic 
for  the  LLBL,  than  the  simple  straight-field  line  model. 


Basic  Equations  and  Model 

The  starting  point  is  the  fully  three-dimensional,  dissipative,  compressible 
magnetohydrodynamic  (MHD)  equations,  written  in  dimensionless  form: 
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where  v(.r,  y,  r,  i )  =  {vr,vy,vz)  is  the  flow  velocity;  B(x,y,z,t)  =  (BT,  By,  Bz)  is  the 
magnetic  field;  R  is  the  Reynolds  number;  Rm  is  the  magnetic  Reynolds  number;  M  is 
the  sonic  Mach  number;  and  Ma  is  the  Alfven  Mach  number.  All  lengths  are  normalized 
by  a  characteristic  perpendicular  half  width  of  the  shear  flow,  a  =  6/2]  the  velocity  by 
the  velocity  jump,  Vo;  the  magnetic  field  by  the  background  magnetic  field,  B0.  The 
sonic  and  Alfven  Mach  numbers  are  defined  as  Ma  =  Vo/v,  and  Ma  =  Vq /?>„,  where 
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vs  =  (7P0/po)l/2  and  va  =  Bo/{4np0)1/2  are  the  sound  speed  and  the  Alfven  speed, 
respectively;  the  Reynolds  number  is  defined  a s  R  =  aVo/v  and  the  magnetic  Reynolds 
number  as  Rm  =  4naVo/{T]c2).  The  pressure,  P,  and  density,  p,  are  normalized  by  their 
values  away  from  the  shear  layer.  Resistive  and  viscous  dissipation  terms  in  the  energy 
equation  are  neglected. 

The  initial  configuration  used  in  this  study  is  a  velocity  shear  layer  in  a  nonuniform 
magnetic  field,  the  nonuniformity  being  created  by  electic  currents  that  connect  the 
plasma  in  the  layer  to  two  conducting  plates  which  serve  to  represent  the  northern 
and  southern  ionospheres.  A  sketch  of  the  configuration  is  shown  in  Figure  1.  The 
quantitative  expressions  describing  the  initial  field  and  plasma  flow  are: 


v(x,y,2,0)  =  tanh(-)  x,  Bx(x,y,z, 0)  =  -a  z  tanh(-), 
la  a 

By(x,  y ,  2,0)  =  0,  Bz(x,  ?/,  2, 0)  =  B0,  p(x,  y,  2, 0)  =  p0, 

.  B2{x,y,z,  0) 

P(x,j/,2,0)  + - — - -  =  const. 


(2) 


The  main  drawback  of  this  configuration  is  that  it  is  not  a  strict  equilibrium.  The  x 
component  of  the  j  x  B  force  is  not  balanced  by  any  other  force  and  therefore  will  result 
in  a  gradual  slowing  down  of  the  flow  and  an  associated  decrease  in  the  velocity  jump 
across  the  layer. 

In  the  present  study,  the  simulation  domain  is  a  rectangular  box  with  periodic 
boundary  conditions  in  the  flow  direction  (x).  The  boundary  conditions  in  the  direction 
across  the  layer  (y)  are  simply  “free-slip”  ones.  We  use  mirror  and  “free-slip”  boundary 
conditions  in  the  equatorial  plane  (2  =  0)  for  all  variables  and  in  the  upper  plane 
for  the  velocity  components.  In  order  to  write  boundary  conditions  for  the  Bx  and 
By  components  of  the  magnetic  field  at  the  upper  or  lower  edge,  we  assume  that  our 
configuration  is  confined  between  two  parallel  resistive  plates  at  2  =  ±LZ.  Using  Ohm’s 
law  and  current  continuity  conditions  at  the  top  of  our  box,  it  is  easy  to  obtain  equations 
for  the  Bx  and  By  components  on  the  upper  plane  at  2  =  Lz.  In  dimensionless  form 
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these  equations  are 


'"•(»■  -  f  S')  -  *4 


dxd2z 
a  dyd2z 


(3) 


where  V^_  =  d2/dx 2  4-  d2/dy 2  and  R^,  is  equal  to  RmyZ/a. 

Finally,  a  boundary  condition  for  the  pressure  can  be  determined  from  the  absence 
of  a  normal  component  of  velocity  on  the  upper  plane. 

For  space  discretization  the  Fourier  pseudo-spectral  representation  was  used  in 
the  x  direction  and  the  Chebyshev  tau  method  in  the  y  and  z  directions.  We  solved 
the  nonlinear  equations  using  an  iteration  scheme  and  a  time  splitting  alternating 
direction  implicit  (ADI)  method,  modified  for  use  with  Fourier  and  Chebyshev  spectral 
discretizations. 


Simulation  Results 

In  order  to  allow  the  KH  instability  to  develop,  a  perturbation  was  imposed  on  the 
initial  configuration  (2).  The  wavelength  of  this  perturbed  mode  is  equal  to  the  entire 
length  of  the  system  in  the  x  direction  and  has  been  chosen  to  be  the  wavelength  of 
the  fastest  growing  mode  (FGM).  Since  we  do  not  intend  to  address  the  question  of 
the  inverse  cascade,  i.e.,  the  formation  of  structures  with  longer  wavelengths  than  those 
predicted  by  linear  theory  [Belmont  and  Chanteur ,  1989],  we  choose  the  system  length 
to  be  equal  to  the  FGM  wavelength:  for  all  our  simulation  runs  the  value  kxa  =  0.45, 
and  hence  Lx  ~  14a,  was  used.  The  amplitude  of  the  initial  perturbation  of  the  y 
component  of  the  velocity  \vy\  was  equal  to  0.01  Vo.  We  also  used  Ly  =  Lz  =  20a, 

R  =  1000,  Rm  =  100,  M  =  0.77  and  Ma  =  1.0  for  all  runs. 

Figure  la  shows  the  temporal  evolution  of  the  maximum  of  the  y  component  of  the 
flow  velocity  |t>y|,  normalized  by  V0,  for  three  different  values  of  namely  0.1,  0.2  and 
0.4,  as  well  as  for  the  purely  two-dimensional  transverse  case  of  the  KH  instability  which 


corresponds  to  R%  =  0.  One  can  see  that  for  all  four  Reynolds  numbers  the  velocity 
perturbation  initially  grows  linearly  and  then  saturates  at  some  level.  The  linear  growth 
rate  decreases  significantly  with  increasing  magnetic  Reynolds  number  Rt.  The  level 
of  the  nonlinear  saturation  is  almost  the  same  for  R%  =0.1  as  for  the  transverse  case, 
/?£  =  0,  but  for  Rz  =  0.4  this  level  is  much  lower. 

Figure  2  shows  streamlines  and  current  distribution,  at  the  end  of  our  simulation, 
for  the  case  of  strong  bending  of  the  magnetic  field  lines  (R%  =  0.4).  The  term  strong 
bending  is  used  only  in  a  relative  sense;  the  maximum  amplitude  of  the  component  of 
the  magnetic  field  parallel  to  the  flow  (x)  is  approximately  20%  of  the  perpendicular  (i) 
component.  The  KH  instability  is  almost  suppressed  in  this  configuration,  the  level  of 
nonlinear  saturation  being  only  2.5  times  higher  than  the  level  of  the  initial  perturbation 
(0.025 Vo).  The  current  sheet  is  only  slightly  disturbed  and  no  localized  current  structures 
are  formed.  We  may  compensate  for  the  decrease  of  the  velocity  jump  across  the  layer 
during  the  simulation,  which  in  this  case  is  rather  big  (AVx0/AVx  =  6),  by  normalizing 
the  amplitude  of  the  perturbed  y  component  of  the  velocity  by  the  magnitude  of  this 
jump.  At  the  end  of  the  simulation,  this  normalized  velocity  component  is  approximately 
0.15,  which  should  be  compared  to  a  value  of  0.33  for  the  transverse  case,  Ry  =  0,  where 
no  decay  of  the  velocity  jump  occurs. 

In  light  of  this  result,  it  seems  likely  that  strong  bending  of  the  magnetic  field 
lines  can  suppress  or  at  least  significantly  slow  down  the  KH  instability.  Therefore, 
previous  two-dimensional  analyses  of  the  stability  of  the  LLBL  [Miura  and  Pritchett, 
1982]  and  MHD  simulations  of  that  layer  [Miura,  1992;  Wei  and  Lee,  1993],  all  of  which 
were  carried  out  in  the  equatorial  plane  and  using  straight  magnetic  field  lines,  may  be 
strictly  valid  only  in  a  limited  region  not  too  far  from  the  subsolar  stagnation  point, 
where  field-line  bending  is  not  strong  enough  to  prevent  the  KH  instability.  But  as  the 
plasma  flow  proceeds  toward  the  tail,  the  curvature  of  the  magnetic  field  lines  increases 
and  one  cannot  neglect  their  tension  any  more. 


Figures  3  and  4  show  current  lines  in  the  system  for  Rz  =  0.2  and  Rz  =  0.1, 
respectively.  The  current  distribution  is  more  complicated  than  in  the  previous  case. 
The  perturbations  introduced  by  the  KH  instability  on  the  initially  uniform  current 
sheet  are  so  strong  that  they  result  in  its  destruction  and  in  the  creation  of  localized 
vortex/current  tubes.  This  process  looks  similar  to  magnetic  island  formation  in 
association  with  the  tearing  mode  instability  of  a  two  dimensional  current  sheet.  Our 
simulations  show  the  formation  of  localized  3D  current  structures  aligned  with  the 
magnetic  field.  However,  the  processes  taking  place  near  the  upper  and  lower  ends 
of  these  tubes  and  the  existence  of  the  inverse  cascade  in  the  system  are  probably 
important  for  understanding  the  long-time  development  of  these  tubes. 

In  order  to  make  geophysical  estimates  from  our  simulations,  the  upper  conducting 
plate  must  be  given  properties  that  mimic  the  ionoshere.  A  relationship  between  the 
conductivity  of  the  plate  and  the  effective  height-integrated  conductivity,  Ept,  of  the 
ionosphere,  based  on  the  conservation  of  magnetic  flux  and  field-aligned  current,  was 
given  by  Sonnerup  [1980].  Using  our  dimensionless  parameter  Rz  =  RmT)T,/a  it  can  be 
written  as: 


c2  1  B  (dxV 
p'  4ir Vo  k  Bi  \dii )  2 


(4) 


The  total  field-aligned  current  in  the  ionosphere,  averaged  over  one  period  in  the  x 
direction,  is  given  by 


<  /i 


_  c  dx  ( 
^  4  x  dii  \ 


(<BX> 

y=0  “  <  Bx  > 

V 

Z=Lz 

V=Ly  I 

z=Lt  / 


(5) 


where  the  subscript  i  denotes  ionospheric  quantities  and  [ dxfdxi ]  is  the  mapping  factor 
for  distances  in  the  main  flow  direction.  Using  typical  values  for  the  magnetospheric 
parameters  [e.g.,  Sonnerup,  ’980],  B  =  3  x  10-8f,  Bi  =  5  x  10-5f,  V0  =  200km/s, 
dx/dx ,  =  51.3  and  taking  the  coupling  factor  k  =  1,  one  can  estimate  the  values  of  Epi 
and  <  /||i  >  as  6.3 Rz  mho  and  1.2 Rz  A/m  respectively.  With  Rz  =  0.1  —  0.4  these 
values  are  roughly  consistent  with  observed  values  of  2-3  mho  and  0.15  A/m. 
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It  should  be  noted  that  many  satellite  observations  have  shown  the  presence  of 
spatially  periodic  bright  spots  in  the  postnoon  auroral  region.  Their  dimensions  usually 
range  approximately  from  50  to  200  km  and  their  separation  is  about  100-500  km 
[e.g.,  Lui  et  ai,  1989].  Recently  Wei  and  Lee  [1993]  suggested  that  these  spots  can 
represent  ionospheric  signatures  of  vortices  created  in  the  LLBL  by  the  KH  instability. 
The  formation  of  vortex/current  tubes  in  our  simulations  indicates  that  regions  of  large 
field-aligned  current  density  and  vorticity,  with  dimensions  comparable  to  observed 
values,  will  occur  in  the  ionosphere  representing  the  footprints  of  these  tubes. 

Conclusion 

The  real  plasma  flow  and  magnetic  field  configuration  in  the  magnetopause¬ 
boundary  layer  region  is  significantly  different  from,  and  far  more  complicated  than,  the 
one  we  use  here.  In  particular,  our  use  of  conducting  plates  to  represent  the  ionosphere 
ignores  Alfven-wave  transit  time  effects.  Nevertheless,  our  model  can  be  used  for  at 
least  a  qualitative  assessment  of  the  role  and  nature  of  the  Kelvin-Helmholtz  instability 
in  the  LLBL  and  other  internal  magnetospheric  shear  layers  with  current  closure  in  the 
ionosphere. 

The  following  main  conclusions  can  be  drawn  from  our  simulations: 

1 .  Magnetic  field-line  bending  leads  to  a  significant  decrease  of  the  growth  rate  of  the 
Kelvin-Helmholtz  instability.  The  value  of  this  decrease  depends  on  the  amount  of 
field-line  bending  and,  hence,  on  the  conductivity  of  the  auroral  ionosphere  to  which 
the  shear  layer  is  magnetically  coupled. 

2.  It  seems  likely  that  sufficiently  strong  bending  of  the  magnetic  field  lines, 
corresponding  in  our  case  to  the  magnetic  Reynolds  number  (based  on  an  equivalent 
ionospheric  conductivity)  /?£  >  0.4,  can  slow  down  the  development  of  the 
Kelvin-Helmholtz  instability  or  suppress  it  altogether. 

3.  A  combination  of  the  tearing  mode  instability  and  the  Kelvin-Helmholtz  instability 
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leads  to  a  formation  of  localized  three-dimensional  vortex/current  tubes.  A 
projection  of  these  structures  into  the  ionosphere  produces  regions  of  enhanced 
field-aligned  current  density  and  vorticity  which  may  represent  auroral  bright  spots 
observed  by  for  example  the  Viking  satellite  [Lui  et  al.,  1989]. 

4.  We  have  also  carried  out  simulations  with  a  velocity  profile  more  representative 
of  the  LLBL,  namely  vr  =  |Vo(l  —  tanh(^)),  the  result  being  an  even  stronger 
suppression  of  the  KH  instability,  presumably  caused  by  the  presence  of  field-line 
curvature  at  the  point  of  inflection  of  the  velocity  profile. 
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Time  evolution  of  amplitudes 


Figure  1.  Temporal  evolution  of  |us|.  Insert:  sketch  of  plasma  flow  and  magnetic  field 
configuration  confined  between  two  conducting  plates 


the  simulation  ( t  =  83a/V0):  (a)  z  —  0.02 Lz\  (b)  z  =  Lz/4 ;  (c)  z  =  Lz!2\  and  (d)  z  =  3Lz/4. 
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Figure  4.  Current  lines  for  Rz  =  0.1,  plotted  at  four  different  times:  (a)  t  =  20a/Vo;  (b) 
t  =  40a/Vo;  (c)  t,  =  60a/Vo;  and  (d)  t  =  82a/V0. 
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An  analysis  is  presented  of  linear  stability  against  tearing  modes  of  a  current 
sheet  formed  between  two  oppositely  magnetized  plasmas  forced  towards  each 
other  in  two-dimensional  steady  stagnation-point  flow.  The  velocity  vector  in 
this  flow  is  confined  to  planes  perpendicular  to  the  reversing  component  of  the 
magnetic  field.  The  unperturbed  state  is  an  exact  resistive  and  viscous 
equilibrium  in  which  the  resistive  diffusion  outwards  from  the  current  sheet  is 
exactly  balanced  by  the  inward  motion  associated  with  the  stagnation-point 
flow.  Thus  the  behaviour  of  the  tearing  mode  can  be  examined  even  when  the 
resistive  diffusion  time  is  comparable  to  or  smaller  than  the  growth  time  of  the 
instability.  The  linear  ordinary  differential  equation  describing  the  mode 
structure  is  integrated  numerically.  For  large  Lundquist  number  S  and  viscous 
Reynolds  number  Re  the  Furth-Killeen-Rosenbluth  scaling  of  the  growth  rate 
is  recovered  with  excellent  accuracy.  The  influence  of  the  stagnation-point  flow 
on  the  tearing  mode  is  as  follows:  (i)  long-wavelength  perturbations  are 
stabilized  so  that  the  unstable  regime  falls  between  a  short-wavelength  and  a 
long- wavelength  marginal  state;  (ii)  for  sufficiently  low  Lundquist  number 
(S  <  12-25)  the  current  sheet  is  completely  stable  to  tearing-mode  perturbations ; 
(iii)  the  presence  of  high  viscosity  reduces  the  growth  rate  of  the  tearing 
instability.  This  effect  is  more  important  at  small  wavelength.  Finally, 
application  of  the  results  from  this  study  to  the  problem  of  solar-wind  plasma 
flow  past  the  earth’s  magnetosphere  is  briefly  discussed. 


1.  Introduction 

Magnetic  reconnection  is  thought  to  be  an  important  process  for  the 
conversion  of  magnetic  field  energy  into  kinetic  and  thermal  energy  in  cosmic 
as  well  as  laboratory  plasmas  (for  reviews  see  e.g.  Vasyliunas  1975;  Sonnerup 
1979  ;  Forbes  &  Priest  1987).  Reconnection  is  initiated  in  thin  current  sheets  as 
a  result  of  the  tearing  mode  (Furth,  Killeen  &  Rosenbluth  1963;  Laval,  Pellat 
&  Vuillemin  1966;  Wesson  1966;  Cross  &  Van  Hoven  1971 ;  Lee  &  Fu  1986),  in 
which  magnetic  islands,  produced  by  non-steady  reconnection,  grow  ultimately 
to  large  amplitudes.  But  current  sheets  observed  in  space  plasmas,  for  example 
in  association  with  solar  flares  or  at  the  earth's  magnetopause,  sometimes 
remain  stable  for  time  periods  far  exceeding  the  growth  time  associated  with 
the  tearing  mode.  Among  the  effects  that  may  influence  tearing  modes  in  a 
*  Now  at  Max-Planck-Inotitut  fur  Extraterrestrische  Phyaik.  8046  Garching,  Germany. 
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substantial  way  are  flows  across  the  current  sheet  (Dobrott,  Prager  &  Taylor 
1977 ;  Killeen  &  Shestakov  1978)  and  flows  along  it  (Bulanov,  Syrovatsky  & 
Sakai  1978;  Einaudi  &  Eubini  1986,  1989;  Chen  &  Morrison  1990;  Ofman,  Chen 
&  Morrison  1991).  However,  a  weakness  of  all  previous  analyses  based  on 
resistive  MHD  is  that  the  unperturbed  state  was  not  an  exact  solution  of  the 
resistive  MHD  equations.  Left  alone,  the  unperturbed  current  sheet  would 
spread  out  with  increasing  time.  Most  of  the  previous  studies  of  the  resistive 
tearing  mode  neglected  this  resistive  diffusion  effect ;  thus  the  results  reported 
are  valid  only  in  the  limit  where  the  diffusion  time  is  much  longer  than  the 
growth  time  of  the  instability.  In  this  limit  the  Lundquist  number  &  =  tD/tA  is 
much  greater  than  unity,  where  tA  is  the  Alfven  transit  time  across  the  current 
layer  and  tD  is  the  resistive  diffusion  time  of  the  current  sheet.  The  Lundquist 
number  of  current  sheets  in  magnetic  reconnection  configurations,  however, 
may  sometimes  be  of  order  unity  (Lee  &  Fu  (1986)  gave  2  <  S  <  100).  The 
studies  of  the  tearing  mode  in  this  regime  have  been  limited  and  not  entirely 
satisfactory :  Lee  &  Fu  (1986)  investigated  the  tearing  mode  in  the  low -<8  regime 
without  taking  into  account  the  resistive  decay  of  the  unperturbed  current 
sheet.  For  a  Harris-type  current  sheet  ( B  oc  tanh  z),  they  found  that  the  growth 
rate  obtained  by  assuming  uniform  conductivity  and  neglecting  the  resistive 
decay  is  essentially  the  same  as  that  obtained  with  a  spatially  varying 
conductivity  of  the  form  a  oc  sech*  x,  which  allows  the  unperturbed  state  to 
remain  time-independent.  However,  this  kind  of  conductivity  profile,  with 
M  ->  0  as  |ar|  -►  oo ,  is  usually  not  relevant  to  the  problem  of  magnetic  reconnection. 
In  order  to  investigate  the  tearing  mode  in  the  low  S- regime  in  a  satisfactory 
way,  one  would  need  to  either  perform  a  stability  .  lalysis  of  a  non-equilibrium 
current  sheet,  taking  into  account  the  resistive  spreading  of  the  sheet  in  a  self- 
consistent  manner,  or  perform  such  an  analysis  on  a  sheet  in  which  the  resistive 
diffusion  is  counterbalanced  by  an  incoming  plasma  flow.  Dobrott  et  al.  (1977) 
recognized  the  importance  of  resistive  decay  of  the  unperturbed  current  sheet. 
However,  although  their  inclusion  of  an  uniform  ‘  diffusion  velocity  ’  to  describe 
the  spreading  of  the  layer  may  give  an  indication  of  the  effect  of  diffusion  on  the 
tearing  mode,  it  is  in  fact  inconsistent  with  the  magnetic  induction  equation.  In 
the  present  paper  we  perform  a  linear  tearing- mode  stability  analysis  of  one 
member  of  the  family  of  exact  resistive  current-sheet  equilibria  found  by 
Sonnerup  &  Priest  (1975)  for  the  case  of  two  oppositely  magnetized  plasmas 
pushed  towards  each  other  in  two-dimensional  stagnation -point  flow.  The 
equilibrium  structure  and  thickness  of  the  unperturbed  current  sheet  is  such 
that  the  resistive  diffusion  outwards  from  the  sheet  is  exactly  balanced  by  the 
inward  motion  associated  with  the  stagnation-point  flow.  The  latter  is  two- 
dimensional,  with  the  flow  confined  to  planes  perpendicular  to  the  reversing 
component  of  the  magnetic  field.  This  equilibrium  allows  us  to  examine  the 
tearing  mode  over  the  entire  range  of  parameter  values. 

The  paper  is  organized  as  follows.  In  §2  the  basic  equations  and  the  relevant 
properties  of  the  equilibrium  configuration  are  reviewed.  In  §3  we  develop  the 
linear  perturbation  equations  and,  by  appropriate  assumptions  concerning  the 
nature  of  the  perturbations,  reduce  them  to  a  form  suitable  for  the  subsequent 
analysis.  In  §4  the  method  of  solution  is  described.  In  §5  numerical  solutions  of 
the  linearized  equations  are  presented.  In  particular,  the  dependence  of  the 
growth  rate  of  the  tearing  mode  on  the  Lundquist  number,  i.e.  the  magnetic 
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Reynolds  number  based  on  the  Alfven  speed,  on  the  viscous  Reynolds  number 
and  on  the  wavelength  of  the  perturbation  is  examined.  Examples  of  the 
eigenmode  structure  are  also  shown.  Finally,  a  discussion  of  the  results  is  given 
in  §6. 


2.  Basic  equations  and  the  equilibrium  state 

The  analysis  is  based  on  the  equations  of  the  incompressible  one-fluid 
resistive  and  viscous  magnetohydrodynamics,  namely 


V.v  =  0, 

(1) 

p^+p(v.Vv)  =  -Vp+i-(VxB)xB  +  *V»v, 

at  fi0 

(2) 

— ^  +  V  x  (v  x  B)  +  ^-^  =  0, 
at  /i0<r 

(3) 

V.B  =  0, 

(4) 

where  B,  v,  p  and  p  are  the  magnetic  field,  plasma  velocity,  density  and  pressure 
respectively,  r/  is  the  dynamic  viscosity  and  <r  is  the  electrical  conductivity, 
both  of  which  are  assumed  uniform  and  constant. 

Steady-state  exact  solutions  to  this  system  of  equations,  for  symmetric  and 
asymmetric  MHD  stagnation-point  flows  in  two  and/or  three  dimensions,  have 
been  given  by  Sonnerup  &  Priest  (1975),  Besser,  Biemat  &  Rijnbeek  (1990)  and 
Phan  &  Sonnerup  (1990).  These  solutions  represent  generalizations  of  the  initial 
work  by  Parker  (1973)  on  resistive  current  layers  in  the  presence  of  two- 
dimensional  stagnation-point  flow.  The  general  form  of  Sonnerup  and  Priest’s 
(1975)  solutions  is 

v0  =  — CiXX+Cjyy  +  CjZZ,  (5) 

B0  =  B9y(x)  y  +£„*(*)  £,  (6) 

where  the  positive  constants  cv  ct  and  c3  are  related  by 

cx  =  Cj  +  C3 

to  ensure  that  the  flow  is  divergence-free.  It  should  be  noted  that  the 
equilibrium  flow,  described  by  (5),  is  irrotational.  Consequently,  the  viscous 
force  on  the  flow  is  zero.  However,  in  the  perturbed  state  the  flow  becomes 
rotational,  and  the  effect  of  viscosity  may  be  important. 

An  early  attempt  to  study  the  behaviour  of  the  tearing  mode  in  general  three- 
dimensional  stagnation-point  flow  and  field  configurations  of  the  type  given  by 
(5)  and  (6)  was  made  by  Sonnerup  &  Sakai  (1981).  They  used  an  analytical 
approach  in  which  the  island  stretching  caused  by  the  accelerating  unperturbed 
plasma  motion  along  the  current  sheet  was  treated  by  use  of  a  method 
introduced  by  Bulanov  et  al.  (1978).  However,  the  analysis  of  Sakai  &  Sonnerup 
remained  unsatisfactory  in  some  aspects  and  was  never  completed.  They  did 
not  consider  the  special  case  c2  =  0,  which  is  much  simpler  than  the  general  case 
on  account  of  the  absence  of  island  stretching.  It  is  this  special  case  that  will  be 
studied  here.  For  c3  =  0  the  flow  is  two-dimensional  and  confined  to  the  (x,  z)- 
plane.  i.e. 

v0  =  c ,  (  xx  -f-  zz),  (7) 
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and  the  magnetic  field  components  B0y  and  B0z  obey 

_L^+e  xdJ*-0 

/i0<r  dx3  1  dx 


P  o 


1  dlB^. 
^~dxr+C'X  dx 


oz  ,  „  dB0t 


+■  Ci  BQt  =  0. 


Also  the  pressure  is  obtained  from  a  Bernoulli-type  equation 

BL+Bl 


2/*0 


(8) 

(9) 

(10) 


where  p0  is  a  constant  of  integration. 

The  resistive  current-sheet  structure  described  by  (8)  and  (9)  has  a 
characteristic  width  of  the  order  of  the  resistive  length 

a  =  (11) 

which  is  obtained  by  equating  the  convective  flow  speed  c,o  to  the  diffusion 
speed  (fi0(ra)~l.  A  basic  property  of  the  current  sheet  obtained  from  (8)  and  (9) 
is  that  the  smaller  the  flow  rate,  i.e.  the  smaller  the  value  of  cx,  the  thicker  will 
be  the  current  sheet.  The  odd  and  even  solutions  of  (8)  and  (9)  can  be  combined 
in  various  ways  in  order  to  produce  the  behaviour  of  B0y  and  B0z  in  the  sheet. 
For  our  purpose  it  will  suffice  to  use  only  the  odd  solution  for  B0y,  which  was 
shown  by  Sonnerup  &  Priest  (1975)  to  be 

B0y  =  AnmxerfKK^i)1*]*  (12) 

where  ±f?max  is  the  magnetic  field  at  x  =  ±  oo.  The  analysis  will  be  valid  for  all 
B0z  satisfying  (9).  The  general  solution  for  B0z  has  been  shown  by  Sonnerup  & 
Priest  (1975)  to  be 


S,=J,|0)W+^|0|/Mr|1  (13) 

where  I(£)  =  exp(-$clp0cr£2). 

The  resulting  equilibrium  plasma  flow  and  B0y,  the  antiparallel  part  of  the 
magnetic  field,  are  shown  in  figure  1.  The  even  and  odd  solutions  of  (13)  are 
shown  in  figure  8.  Note  that  these  B0z  solutions  vanish  as  |x|->  oo. 
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3.  Linear  perturbations 

Returning  to  (l)-(3),  we  now  introduce  small  perturbations  v1  and  Bt  in  the 
velocity  and  magnetic  field  so  that 

v  =  v0+v,(x,  y,  I), 

B  =  B0-f  B^x,  y,<)- 

Note  that  these  perturbations  are  taken  to  be  independent  of  the  co-ordinate  z. 
Thus  the  wave  vector  of  the  tearing  mode  is  assumed  to  be  along  the  y  axis. 
Equation  (1)  is  satisfied  identically  by  writing 

Vj  =  V^xz  +  zr  u(x,y,t),  (14) 

where  is  the  stream  function  for  the  x  and  y  components  of  v,. 

The  pressure  may  be  eliminated  from  the  problem  by  taking  the  curl  of  (2), 
the  result  being 

V  x  (v  x  fl)  =  —  V  x  [(V  x  B)  x  B]  +  ^V*a,  (15) 

dt  n9p  p 

where  12  is  the  vorticity,  i.e. 

flsVxv  =  Vxv!  =  Vt>lfxz— zVV-  (16) 

The  last  two  equalities  in  (16)  follow  from  the  fact  that  V  x  v0  =  0  and  from  (14) 
respectively. 

We  now  examine  the  z  component  of  the  vorticity  equation  (15),  and  the  x 
component  of  the  induction  equation  (3).  The  linearized  versions  of  these 
equations  are 


=  *S*(V* 

N  P\ 


1  d*B, 


dx* 


'l* 


(17) 


and 


(18) 


respectively.  It  should  be  noted  that  the  only  perturbation  quantities  contained 
in  the  above  equations  are  ijr  and  Blz.  Thus  we  see  that  it  will  suffice  to  analyse 
these  two  equations  to  determine  the  evolution  of  the  perturbations.  Once  \fr 
and  Blx  are  known,  the  perturbed  quantities  vu  and  Bu  may  subsequently  be 
obtained  from  the  z  components  of  (2)  and  (3),  namely 


“ci 


a  d  1 

a.  Cj  C|X  - 

dt  ox  pt<r 


u  HoP\  °*  dy 

♦*.*) 

(19) 

—B  ^ t>u 

V  )0U  dy 

5xlfdBu 
dy  dx 

(20) 

It  should  be  noted  that  B0t  does  not  enter  into  (17)  and  (18) ;  thus  it  does  not 
play  a  role  in  driving  the  instability.  It  does,  however,  entrain  vu  and  Bu  via 
( 19)  and  (20).  For  B0z  =  0  the  appropriate  solutions  of  these  latter  equations  are 
vu  =  0  and  Bu  =  0. 
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Solutions  to  (17)-(20)  are  now  sought  in  the  form 


*  1 

f  #(X)  ' 

*«. 

Blx(x) 

vu 

>Bu(z), 

exp(iky+yt), 


(21) 


where  lc  and  y  are  the  wavenumber  and  growth  rate  respectively.  Substitution 
of  these  expressions  into  (17)— (20)  gives 


| y+Cj-e,; 


dx  n9p' 


d 

1  / d * 

dx 

ti0<r\dz* 

(22) 

(23) 

•jbi***-**"^)- 

(24) 

(25) 

We  now  non-dimensionalize  (22)— (25)  by  introducing  the  following  dim¬ 
ensionless  variables: 


**  =  -,  k*  —  ka,  y*  — 


B, 


»  -  rsa 

Here  a  is  the  resistive  length  given  by  (11),  vA  is  the  Alfven  speed  based  on  the 
maximum  magnetic  field  Bm„,  and  c,o  is  the  characteristic  diffusion  speed. 
Thus  the  quantity  S,  known  as  the  Lundquist  number,  is  the  magnetic 
Reynolds  number  cravA  based  on  the  Alfven  speed  and  the  diffusion  length. 
It  should  also  be  noted  that  S  =  where  M A  is  the  Alfven  Mach  number 
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at  the  edges  of  the  current  sheet  (i.e.  at  x  =  ±a).  The  quantity  Re  is  the  viscous 
Reynolds  number  based  on  the  characteristic  diffusion  speed  and  the  diffusion 
length.  Finally,  the  renormalization  of  Bu  is  such  that  it  is  phase-shifted  by  90° 
relative  to  the  other  perturbed  quantities. 

Substitution  of  the  above  expressions  into  (22)— (25)  gives 


,  m_d  1  /  d' 
X  dx*  Re\dx*x 


(26> 

{?’+'-z’£i-£r>+k”)B"  =  (27) 

[r*+>— ,28, 


while  the  equilibrium  magnetic  field  profiles,  given  by  (12)  and  (13),  become 

B*y  =  erf  (2-***),  (30) 

&tz  =  £?*(<>)  exp  ( -  Jx**) + (0)  exp  ( -  £x**)  £  dg exp  (|£*).  (31 ) 

The  coupled  equations  (26)-(29)  describe  all  small-amplitude  MHD  wave  modes 
of  the  current  sheet  that  have  their  propagation  vector  along  the  y  direction. 
We  shall  examine  only  the  tearing  mode,  for  which  the  amplitudes  B*z(x*)  and 
\Jf*(x*)  are  respectively  even  and  odd  functions  of  x*,  and  for  which  these 
quantities  decay  rapidly  with  increasing  distance  |x*|  from  the  centre  of  the 
current  sheet.  The  parities  of  vu  and  Bu,  on  the  other  hand,  depend  on  B0t.  In 
particular,  Bu  has  the  same  parity  as  B0t,  whereas  vu  has  the  opposite  parity. 
The  perturbed  quantities  vu  and  Bu  are  also  required  to  vanish  at  large 
distance.  The  procedure  for  solving  (26)— (29)  is  as  follows.  Equations  (26)  and 
(27)  are  not  coupled  to  the  other  equations,  and  may  be  solved  for  y*,  B*x(x*) 
and  y*(x*)  first.  Solutions  for  vf,(x*)  and  B*z{x*)  are  subsequently  obtained 
from  (28)  and  (29). 

In  this  paper  we  determine  the  largest  real  growth  rates  of  the  tearing  mode 
by  a  trial-and-error  method  of  solution.  In  the  case  analysed  by  Furth  et  al. 
(1963)  it  was  determined  that  no  overstable  tearing  modes  occur.  Although  the 
eigenvalue  we  find  is  purely  real  and  agrees  with  the  results  of  Furth  et  al.  for 
large  S ,  we  have  not  found  a  way  to  prove  rigorously  that  the  eigenvalue  with 
the  largest  real  part  is  purely  real.  The  existence  of  overstable  modes  therefore 
remains  an  open  question  in  the  case  investigated  here. 


4.  Method  of  solution 

In  contrast  with  the  standard  boundary -layer  approach  (Furth  et  al.  1963), 
which  involves  matching  of  solutions  from  an  outer  and  an  inner  region,  we 
solve  (26)-(29)  directly  over  the  entire  region.  This  approach  allows  us  to 
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investigate  the  entire  range  of  Lundquist  number  S,  in  particular  values  of  S  of 
order  unity.  Since  only  the  odd  solutions  for  ft*(x*)  and  the  even  solutions  for 
5fx(jr*)  are  sought,  the  equations  need  only  be  solved  over  half  of  the  range  of 
x*.  for  example  for  x*  ^  0.  The  following  boundary  conditions  are  imposed 
when  solving  (26)  and  (27): 

ft*  =  ft*"  =  0,  B*x  =1  at  x*  =  0,  j 

ft*  =  ft*"  =  B*z  -*■  0  as  x*  oo .  J 

The  growth  rate  y*  is  then  determined  by  requiring  the  solution  for  B*x(x*) 
to  be  an  even  function  of  xm  ;  that  is,  B*'(x*  =  0)  =  0.  The  method  for  obtaining 
the  eigenvalue  y*  is  as  follows :  we  start  with  a  trial  value  of  y* ;  the  value  of 
Bf} '(x*  —  0)  is  then  obtained  by  solving  (26)  and  (27),  In  general,  B*x(x*  =  0) 
does  not  vanish.  The  value  of  y*  is  subsequently  adjusted  until  B*x(x *  =  0) 
becomes  zero.  The  resulting  y*,  ft*  and  B*x  are  then  the  eigenvalue  and 
eigenfunctions  of  the  tearing  mode.  The  method  is  somewhat  similar  to  that 
used  by  Wesson  (1966). 

Once  y*,  B*x(x*)  and  ft*(x*)  have  been  determined,  one  can  then  proceed  to 
obtain  solutions  for  vft(x*)  and  B*t(x*)  from  (28)  and  (29),  for  a  given  B%z(x*) 
satisfying  (9),  subject  to  the  condition  that  these  perturbed  quantities  vanish 
as  |a;*|  ->  oo.  The  magnitudes  of  vft(x*)  and  B*t(x*)  are  then  determined  uniquely 
by  the  solutions  for  ft*,  B*x,  y*  and  B$z. 

The  coupled  ordinary  differential  equations  (26)-(29)  are  solved  using  a 
finite-difference  method  with  variable  mesh  size.  A  large  number  of  com¬ 
putational  meshes  are  concentrated  near  x*  —  0,  where  the  solutions  of  ft*;  B*x, 
v*z  and  B*z  display  their  steepest  gradients.  In  the  actual  calculation  the  outer 
boundary  x^  is  located  at  large  but  finite  distance  from  the  origin.  The  choice 
of  xx  is  such  that  its  location  does  not  affect  the  resulting  y*. 


5.  Results 

The  growth  rate  and  the  eigenmode  structure  described  by  ft*(x*),  B*x(x*), 
vfz(x*)  and  Bfz(x*)  depend  on  the  parameters  S,  Re  and  k*.  Figure  2  displays  the 
dependence  of  the  normalized  growth  rate  y*  on  the  normalized  wavenumber 
k*  for  several  values  of  the  magnetic  Reynolds  number  S  and  for  large  viscous 
Reynolds  number  (Re  —  10*).  For  each  5  there  is  a  maximum  growth  rate  y*., 
and  a  corresponding  wavenumber  k^^.  This  maximum  growth  rate  occurs  at 
longer  wavelength  (smaller  wavenumber)  for  larger  S.  For  S  <  &crltlca,  %  12-25 
the  maximum  growth  rates  are  negative,  i.e.  all  perturbations  are  damped.  For 
S  >  <ScrltlC4l  the  unstable  regime  falls  between  a  large- wavenumber  ^Jpper  (short- 
wavelength)  marginal  state  (y*  =  0)  and  a  small-wavenumber  &*ower  (long- 
wavelength)  marginal  state.  The  dependence  of  k*vptT,  fc,*wer,  and  on  S  is 
shown  in  figure  3.  The  existence  of  a  long-wavelength  marginal  state  may  also 
be  deduced  by  examining  the  governing  equations  in  the  limit  k*  =  0.  In  this 
limit  (27)  becomes 


d*B*  dB* 

- 4.  x* - 

dx*%  dx* 


(y*+L)B*  =  0. 


(33) 


This  equation  admits  solutions  that  are  even  functions  of  x*  and  that  vanish  as 
|x*|  oo  if  and  only  if  y*+ 1  is  negative,  i.e.  only  if  y*  <  —  1 .  Thus  the 
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Figure  2.  Normalized  growth  rate  y*  =  y/cl  as  a  function  of  normalized  wavenumber  k*  — 

lea  for  Lundquist  numbers  S  *  5  ( A — A),  7-5  (□ — □),  10(9 — 9).  *2-25  ( - ),  15(0 — O). 

20  (O — O).  50  (■ — ■)  and  100  (  x  —  x  )  and  for  viscous  Reynolds  number  Re  *  10*.  The 
growth  rates  are  negative  for  S  <  12*25. 

configuration  is  stable  for  k*  —  0.  Since  (for  sufficiently  large  S)  it  is  unstable  for 
some  k *,  a  marginal  state  must  exist  at  long,  but  finite  wavelength. 

Figure  3  shows  again  that  the  unstable  regime  narrows  as  S  decreases  and 
converges  to  a  point  as  S  reaches  ScrlUcal.  Below  Scrixlc%1  the  current  sheet  is 
stable  to  tearing-mode  perturbations.  For  large  S,  k*pper  and  k*ower  asympto¬ 
tically  approach  0*733  and  zero  respectively,  while  obeys  approxi¬ 
mately  the  following  relation : 

*m»x  *  0*91  S-**45  (104  <  S  <  5  x  10*).  (34) 

It  should  be  pointed  out  that  the  asymptotic  value  of  k*pper  may  be  obtained 
by  performing  standard  boundary-layer  analysis  of  the  type  first  performed  by 
Furth  et  al.  (1963)  for  a  hyperbolic  tangent  magnetic  field  profile.  In  that 
analysis  the  short-wavelength  marginal  state  is  found  by  setting  pper) 

=  0,  where  AoUter  is  the  jump  in  B*~l  dB*/dx *  of  the  outer  solution  across  the 
singular  layer.  Our  boundary-layer  analysis  for  an  error-function  magnetic  field 
profile  gives  k*pp^r  =  0*733. 

Figures  4 (a,  6)  show  the  maximum  growth  rate  as  a  function  of  S  for 
S  >  100  and  S  $  100  respectively.  For  large  S  the  curve  is  almost  a  straight  line 
on  a  log-log  scale,  and  may  be  approximated  by 

y£«  *  0-39S*608  ( 104  <  S  <  5  x  10*).  (35) 
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Figure  3.  Wavenumbers  of  fastest  growth,  1^,,  of  short- wavelength  marginal  state 
an(*  ^iow«r  °f  long-wavelength  marginal  state  as  functions  of  Lundquist  number  S  =  fi0  <ravM 
for  viscous  Reynolds  number  Re  -  10*.  The  unstable  k*  regime  narrows  as  S  decreases. 
Below  <ScrlUc%1 «  12-25  the  current  sheet  is  stable  to  tearing-mode  perturbations.  For  S  >  104 
the  curve  may  be  approximated  by  as  0-91/S"**4*. 


Figure  4.  Maximum  growth  rate  yZ..  as  a  function  of  Lundquist  number  S  for  Re  =  10*  and 
S  >  100  (a)  and  S  ^  100  (6).  For  S  >  10*  the  curve  may  be  approximated  by  y*„  «  0-39S****. 
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Figure  5.  Normalized  growth  rate  y *  **  y/c,  as  a  function  of  normalized  wavenumber  k*  — 
ka  for  viscous  Reynolds  numbers  Re  =  1  (□ — □),  10  (ff — M)>  100  (O — O).  10* 
(O — O)  and  104  ( x —  x  )  and  for  Lundquist  number  S  *  50. 

The  exponents  in  (34)  and  (35)  are  very  close  to  the  asymptotic  values  (as 
S-*  oo)  of  £  and  J  respectively,  obtained  by  Furth  et  al.  (1963)  in  their  study  of 
the  tearing  mode  without  equilibrium  flow,  indicating  that  the  effect  of  the 
equilibrium  flow,  other  than  that  of  determining  the  current-sheet  width,  is 
negligible  when  S  is  large.  For  small  S,  however,  the  relation  between  and 
S  differs  considerably  from  (35),  and  also  does  not  agree  with  Lee  &  Fu’s  (1986) 
result  for  low  S.  In  particular,  those  authors  did  not  find  a  stable  regime  for  low 

5. 

The  dependence  of  the  growth  rate  on  viscosity  is  illustrated  in  figures  5  and 

6.  Figure  5  shows  the  normalized  growth  rate  y*  as  a  function  of  the  viscous 
Reynolds  number  Re  for  5  =  50.  It  is  seen  that  the  viscosity  has  a  stabilizing 
effect.  This  effect  is  more  important  at  short  wavelength  (large  k*)  and  when  the 
viscous  Reynolds  number  Re  is  small  (Re  <  100).  It  should  also  be  noted  that 
the  fastest-growing  wavelength  increases  with  increasing  viscosity  (decreasing 
Re).  The  stabilizing  effect  of  the  viscosity  may  also  be  seen  in  figure  6,  where 
‘-’critical  »s  shown  as  a  function  of  Re.  It  should  be  noted  that  ScrlUcft,  decreases 
with  increasing  Re.  For  Re  >  103,  *ScrlUcsl  is  close  to  its  asymptotic  value  of 
12-25. 

In  figures  7  (a,  b)  we  show  two  examples  of  the  eigenmode  structure  in  terms 
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Figvre  7.  Equilibrium  magnetic  field  B^(x*)  and  perturbed  eigenmodes  <lr*(x*)  and  Bfjx*) 
as  functions  of  x*  —  x/a  for  Re  =  10*.  k*  =»  04,  and  (a)  5  =  5  and  (6)  S  =  1000. 

of  \}f*(x*)  and  B*x{x*)  for  S  =  o  and  1000  respectively,  and  for  Re  =  10s  and 
fc*  =  0*4,  illustrating  the  difference  in  the  behaviour  of  the  eigenmodes 
depending  on  the  S  regime.  In  the  S  =  1000  case  the  shapes  of  the  eigenmodes 
are  similar  to  those  obtained  in  previous  studies  of  the  tearing  mode  with  high 
S,  and  where  the  effect  of  resistive  decay  of  the  current  sheet  is  neglected  (see 
e.g.  Killeen  &  Shestakov  1978).  In  particular,  an  ‘inner  region ’  appears  that  is 
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Figure  8.  Equilibrium  magnetic  field  ££(z*)  and  perturbed  eigenmodes  vf,(x*)  and  B%t(x*) 
as  functions  of  x*  =  x/a  for  even  (a,  b)  and  odd  (c,  d )  magnetic  field  for  S  =  5  (a,  c)  and 
S  —  1000  (6,  d).  The  viscous  Reynolds  number  and  wavenumber  are  Be  —  10*  and  k*  «  0-4. 

much  thinner  than  the  current-sheet  width.  In  the  S  =  5  case,  as  expected,  the 
‘  inner  region  ’  is  comparable  to  the  resistive  length.  In  this  regime  the  results 
of  the  previous  studies,  which  excluded  the  diffusion  effect  of  the  equilibrium 
magnetic  field,  are  not  valid.  Note  in  particular  that  the  tearing  mode  is 
damped  for  S  =  5  even  though  the  curve  for  B*x(x*)  displays  the  dimple  at 
x*  —  0  normally  associated  with  unstable  behaviour.  For  convenient  com¬ 
parison  of  length  scales  the  equilibrium  magnetic  field  B*y  is  also  shown  in 
figure  7. 

Figures  8  (a-d )  display  the  eigenmode  structures  in  vft(x*)  and  B*t(x*)  for  the 
same  cases  as  in  figure  7,  and  for  a  purely  even  (figures  8a,  6)  and  a  purely  odd 
(figures  8c,  d )  equilibrium  magnetic  field  B0z.  For  an  arbitrary  satisfying  (9) 
the  resulting  eigenmode  structures  in  v*t(x*)  and  B*,(x*}  are  linear  combinations 
of  the  even  and  odd  eigenfunctions  shown.  Note  that  the  even  solutions  of 
display  large  local  curvature,  for  example  x*  —  0  in  figure  8(c).  This  behaviour 
is  a  direct  consequence  of  the  large  Re  value.  Note  also  that  the  spatial  scales 
are  different  in  figures  7  and  8. 
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6.  Discussion  and  conclusion 

We  have  examined  the  linear  stability  of  a  current  sheet  against  the  tearing 
mode  in  the  presence  of  equilibrium  stagnation- point  flow,  the  latter  being 
confined  to  the  plane  perpendicular  to  the  lc  vector  of  the  tearing  mode.  The 
entire  range  of  Lundquist  numbers  S  and  viscous  Reynolds  numbers  Re,  where 
S  =  fiQ<ravA  and  Re  =  claip/tf,  has  been  explored.  Our  main  findings  can  be 
summarized  as  follows. 

(i)  Long-wavelength  perturbations  are  stabilized  by  the  stagnation-point 
flow  so  that  the  unstable  regime,  if  it  exists,  falls  between  a  short- wavelength 
and  a  long-wavelength  marginal  state. 

(ii)  For  large  S  (>  104)  the  Furth  et  al.  (1963)  scaling  of  the  growth  rate, 
y*,,  oc  and  of  the  wavenumber  at  maximum  growth  rate,  fc*,w  oc  S~K  are 
approximately  recovered,  indicating  that  the  effect  of  the  equilibrium  flow, 
other  than  that  of  determining  the  current-sheet  width,  is  small  when  S  is  large. 
In  dimensional  form  the  growth  rate  and  the  wavenumber  can  be  expressed  as 
ymM  oc  i^A  c\  erf  and  kmMX  oc  <A  respectively.  These  relationships  show  that  a 
faster  flow  (larger  ct)  and/or  a  larger  conductivity  <r  reduce  the  current-sheet 
width,  and  thereby  enhance  the  growth  rate  of  the  instability,  at  the  same  time 
decreasing  the  wavelength  at  which  maximum  growth  occurs. 

(iii)  For  small  S  the  stabilizing  effect  of  the  equilibrium  flow  is  evident.  When 
S  <  <ScrlUcml  as  12-25,  the  current  sheet  is  completely  stable  to  all  tearing-mode 
perturbations,  a  result  not  obtained  in  the  absence  of  the  equilibrium  flow. 

(iv)  The  viscosity  has  the  effect  of  reducing  the  growth  rate.  This  effect  is 
more  noticeable  at  short  wavelengths  and  in  flows  with  small  viscous  Reynolds 
number  Re.  As  a  result,  £crlUcal  increases  with  smaller  Re. 

(v)  The  mode  structure  in  the  (x,y)  plane  and  the  growth  rate  of  the 
instability  remain  unaffected  by  the  presence  of  a  non-zero  magnetic  field 
component  B0z(x)  satisfying  (9).  But  when  B^  4=  0  the  mode  structure  includes 
perturbations  in  the  field  and  flow  components  in  the  z  direction. 

The  main  difference  between  our  study  and  previous  investigations  of  the 
resistive  tearing  mode  is  that  we  have  started  from  an  exact  equilibrium 
current-sheet  configuration  in  which  the  resistive  widening  of  the  sheet  is 
exactly  counterbalanced  by  the  stagnation-point  flow.  It  is  this  feature  of  the 
unperturbed  current  layer  that  allows  us  to  investigate  the  entire  range  of 
Lundquist  numbers  S,  in  particular  values  of  S  of  order  unity.  Our  results  for 
low  S  differ  significantly  from  those  obtained  by  Lee  &  Fu  (1986),  who  did  not 
use  an  exact  unperturbed  equilibrium  in  their  analysis. 

It  may  be  thought  that  the  flow  geometry  studied  in  this  paper  has  limited 
practical  applications.  However,  it  has  been  argued  (Pudovkin  &  Semenov 
1977  a,  6;  Sonnerup  1980)  that,  in  the  absence  of  reconnection,  steady-state  flow 
of  a  magnetized  highly  conducting  plasma  past  a  diamagnetic  object  will  lead 
to  the  formation  of  a  stagnation  line  rather  than  a  stagnation  point  on  the 
upstream  face  of  the  object.  This  stagnation  line  is  aligned  with  the  magnetic 
field  embedded  in  the  impinging  plasma  flow,  so  that  the  unperturbed  flow 
configuration  becomes  similar  to  that  examined  here.  One  application  may  be 
the  flow  of  solar-wind  plasma  past  the  earth’s  magnetosphere,  where  a 
stagnation  line  may  be  formed  near  the  subsolar  point  of  the  magnetopause.  To 
apply  our  results  to  the  magnetopause,  we  must  estimate  the  Lundquist 
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number  based  on  the  thickness  of  the  magnetopause  current  layer.  If  we  adopt 
the  value  of  (/i0<r)-1  =  10*  s  m-8  (see  e.g.  Sckopke  et  al.  1981)  as  an  upper  limit 
for  the  resistive  diffusion  coefficient,  and  use  lower  limits  of  100  km  for  the 
thickness  of  the  current  layer  and  100  km  s-1  for  the  Alfven  speed  just  outside 
the  current  layer,  we  obtain  a  lower  limit  of  S  =  <ravA  %  10.  This  value  is 
consistent  with  Lee  &  Fu’s  (1986)  estimate  of  2  <  S  <  100.  The  fact  that  in  our 
study  the  current  sheet  is  found  to  be  stable  for  S  <  12-25  therefore  suggests 
that  the  magnetopause  current  layer  may  at  times  be  stable  or  only  weakly 
unstable  to  tearing-mode  perturbations. 

This  research  was  supported  by  the  National  Science  Foundation,  Atmo¬ 
spheric  Sciences  Division  under  Grant  ATM-8807645  and  by  the  Air  Force 
Geophysics  Laboratory  under  Contract  F19628-90-K-009  to  Dartmouth 
College.  The  study  was  inspired  by  research  into  the  problem  of  tearing  modes 
in  the  presence  of  three-dimensional  stagnation-point  flow  performed  by 
Professor  J.-I.  Sakai  during  an  extended  stay  at  Dartmouth  College  in 
1980-1981.  The  authors  are  grateful  to  Dr  R.  Richard,  L.  N.  Hau  and  W. 
Lotko  for  helpful  discussions. 
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MAGNETIC  HELD  MAXIMA  IN  THE  LOW  LATITUDE  BOUNDARY  LAYER 
B.  Sonnemp'  ,  G.  Paschmann-  ,  T.-D.  Phan-,  and  H.  Liihr^ 


Abstract.  The  magnetic  field  often  exhibits  a  maximum  in 
the  Earth's  low-latitude  boundary  layer.  We  show  examples 
of  this  behavior,  using  data  from  the  AMPTE/1RM  spacecraft, 
and  argue  that  two  fundamentally  distinct  causes  exist  for  the 
excess  field:  (i)  a  depression,  within  the  layer,  of  the 
population  of  medium-energy  ions  of  magnetospheric  origin; 
(ii)  field  curvature  effects  associated  with  undulations  of  the 
magnetopause  itself. 

1.  Introduction 

The  frequent  presence  of  a  magnetic-field  strength 
maximum  near  the  magnetospheric  edge  of  the  magnetopause 
(MP)  has  been  noted  by  Neugebaueret  al.  [1974],  Sonnerup 
and  Ledley  [1979],  and  Hall  ct  al.  [1991].  Sometimes  the 
region  of  field  enhancement  has  an  inner  edge  that  coincides 
with  the  eanhward  edge  of  the  low-latitude  boundary  layer 
(LLBL)  and  beyond  which  the  field  magnitude  stays  nearly 
constant  at  its  magnetospheric  level;  sometimes  the  field 
decays  from  its  maximum  to  its  magnetospheric  level  in  a  more 
gradual  manner,  as  one  moves  from  the  inner  edge  of  the 
LLBL  into  the  magnetosphere  proper.  We  shall  refer  to  these 
two  pure  classes  of  field  behavior  as  Type  I  and  Type  II, 
respectively,  although  we  stress  that  they  represent  an 
oversimplification:  more  commonly,  a  mixture  of  the  two  is 
seen.  At  first,  the  occurrence  of  field  maxima  seems 
paradoxical.  After  all,  the  region  just  earthward  of  the  MP  is 
usually  occupied  by  the  LLBL,  a  narrow  region  of  dense 
magnetosheath-like  plasma  flowing  along  the  MP  more  or  less 
in  the  antisolar  direction.  One  would  expect  the  diamagnetic 
effect  of  this  plasma  to  produce  a  field  depression  rather  than  a 
field  enhancement  within  the  LLBL,  a  feature  that  is  in  fact 
also  seen  when  the  LLBL  density  is  high. 

Neugebauer  et  al.  [  1974]  did  not  discuss  the  LLBL  but  they 
suggested  that  the  field  maxima  may  be  caused  by  loss  of 
energetic  magnetospheric  particles,  whose  gyromotion  carry 
them  up  to  the  MP  where  they  escape  to  the  magnetosheath. 
However,  they  lacked  the  instrumentation  needed  to  check  this 
hypothesis.  Here,  we  use  data  from  the  AMPTE/IRM 
spacecraft  to  demonstrate  that  field  enhancements  of  Type  I  are 
in  fact  colocated  with  the  LLBL  and  that  the  densities  of 
medium-energy  magnetospheric  ions  (9<E<40  keV)  and 
electrons  (2<E<40  keV)  are  indeed  depressed  in  them.  Thus 
there  is  partial  agreement  with  the  Neugebauer  et  al. 
explanation:  the  excess  magnetic  pressure  compensates  for  the 
defect  in  plasma  pressure  caused  by  the  absence  of  energetic 
magnetospheric  particles.  However,  in  that  explanation, 
electrons  with  their  much  smaller  gyroradii  should  drop  out 
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very  close  to  the  MP  rather  than  eanhward  of  the  inner  edge  of 
the  LLBL,  as  we  observe  them  to  do.  In  Type  II  events,  there 
is  a  gradual  loss  of  high  energy  panicles  (40<E<400  keV)  as 
the  MP  is  approached  but  no  significant  loss  of  medium- 
energy  parades  in  most  of  the  field-enhancement  region. 

Using  AMPTE/UKS  data.  Hall  et  al.  ]1991]  found  a 
minimum  in  total  electron  pressure  precisely  where  the 
magnetic  field  overshoot  occurs:  they  referred  to  this  region  as 
a  “depletion  layer”  (we  reserve  this  term  for  the  region 
immediately  outside  the  MP  where  plasma  depletion  and  an 
associated  field  overshoot  is  sometimes  seen).  This  minimum 
in  electron  pressure  is  a  direct  result  of  the  coexistence  of  two 
electron  populations  in  the  LLBL:  cool  magnetosheath 
electrons  decreasing  in  density  with  increasing  distance  inward 
from  the  MP;  and  hot  (but  far  more  tenuous)  magnetospheric 
electrons,  decreasing  in  density  with  increasing  distance 
outward  from  the  inner  edge  of  the  LLBL.  We  emphasize  that 
the  electrons  usually  only  play  a  minor  role  in  the  overall 
pressure  balance  across  the  MP/LLBL:  the  major 
contributions  to  Ptot  =  (P±  +  B^!2\i0)  come  from 
magnetosheath  ions,  medium-enegy  magnetospheric  ions  and 
the  magnetic  field.  The  imports  inference  between  Type  I 
and  Type  II  events  is  that  Ptot  is  approximately  constant  from 
the  magnetosheath  through  the  MP/LLBL  into  the 
magnetosphere  in  Type  I  whereas  a  sometimes  large  maximum 
in  Ptot  occurs  earthward  of  the  MP  in  Type  II.  This  maximum 
is  caused  by  excess  magnetic  pressure. 

2.  Type  I  Event 

Figure  1  shows  a  Type  I  crossing,  i.e.,  a  case  where  the 
magnetic  field  overshoot  is  confined  to  the  LLBL.  This 
inbound  pass  through  the  MP/LLBL  region  on  July  3,  1985. 
took  place  on  the  dusk  flank  of  the  magnetosphere  ( 1930h  LT) 
near  the  equatorial  plane  (-6.9°  GSE  latitude)  at  a  geocentric 
distance  of  17  Re*  The  motion  of  the  MP  relative  to  the 
spacecraft  was  complicated.  We  go  through  the  event 
backwards,  i.e.,  starting  at  the  right-hand  edge  of  the  figure,  at 
1 625  UT,  where  the  spacecraft  was  in  the  magnetosphere,  in  a 
more  or  less  uniform  field  of  B  =  15  nT  (Pb  £  0.09  nPa. 
panel  6),  with  total  ion  density,  Np,  {E<40  keV,  which  is 
mainly  low  energy  ions)  as  well  as  densities  of  medium- 
energy  ions,  N2p  (9<E<40  keV),  and  electrons,  N 2e 
(2<E<40  keV),  (panel  1)  at  their  magnetospheric  levels.  The 
bulk  speed,  Vp,  of  the  magnetospheric  plasma  was  modest 
(panel  2),  its  temperatures,  Tp  and  Te, were  high  (panel  3). 
Going  backwards  in  time,  the  first  indication  of  the  LLBL 
being  approached  occurred  around  1621:50  UT  where  N2e 
and  Te  started  to  drop.  At  1621:10  UT,  the  LLBL  proper  was 
entered:  Np  increased  abruptly,  while  Tp  and  Te  dropped 
equally  abruptly  to  intermediate  levels  characterizing  the 
LLBL.  The  medium-energy  ion  density,  N2p ,  dropped  more 
gradually.  The  field  magnitude.  B,  rose,  also  gradually  at 
first,  and  then  more  rapidly  until  a  plateau  was  reached  at  B  = 
25  nT  (Pb  =  0.25  nPa,  panel  6).  The  MP  encounter  occurred 
at  1619:55  UT.  Here  B  dropped  abruptly  to  less  than  5  nT 
(PB  <  0.01  nPa);  simultaneously,  the  field  direction  changed 
(panels  4  and  5),  Tp  and  Te  decreased  while  Vp  and  Np  both 
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Fig.  1.  AMPTE/IRM  data  during  inbound  pass  through  the 
MP/LLBL  region  on  July  3,  1985.  Panel  1:  proton  number 
density,  Np  (cm'3;  E<40  keV);  medium-energy  proton  and 
electron  number  densities,  N2p  (cm- 3;  9<E<40  keV)  and 
/V2e(cm'3;  2<E<40  keV).  Panel  2:  plasma  bulk  flow  speed, 
Vp  (km/s).  Panel  3:  proton  and  electron  temperatures,  Tp  and 
Te.  (106°K;  E<40K  keV).  Panels  4  and  5:  magnetic  field 
azimuth  angle,  cp£,  and  elevation  angle,  Xfl,  (degrees). 
Angles  refer  to  boundary-normal  coordinates,  LMN,  with  X  = 
0  in  the  MP  tangent  plane  and  X>0  for  an  outward  field;  also, 
<pfl  =  0  along  +L  and  tp  g  =  90°  along  +M.  Panel  6:  total 
perpendicular  pressure  Ptot  =  Pp  +  Pc  +  Pb  and  magnetic 
pressure  Pp  =  B-! 2p0  (nPa).  Panel  7:  proton  and  electron 
perpendicular  pressures.  Pp  and  Pc.  (nPa;  E<40  keV). 


increased  abruptly  to  their  values  in  the  magnetosheath. 
Continuing  backward  in  time,  several  additional  full  or  partial 
MP  crossings  occurred:  we  do  not  discuss  these  in  detail, 
except  to  note  the  similarity  between  the  MP/LLBL  encounter 
just  discussed  and  that  between  1600:50  and  1603:00  UT. 

The  most  remarkable  feature  in  Figure  1  is  the  approximate 
constancy,  over  a  30  minute  period,  of  the  total  pressure, 
P tot-  (panel  6,  upper  curve),  in  spite  of  large  temporal 
variations  in  the  magnetic  pressure  (panel  6,  lower  curve). 
The  large  ion  (E<40  keV)  and  electron  (E<40  keV) 
perpendicular  pressure  variations.  Pp  and  Pg,  that  compensate 
for  tne  changes  in  magnetic  pressure  are  shown  explicitly  in 
pa.'.el  /.  Note  the  deep  minimum  in  Pc  during  the  LLBL 


encounters  (at  -1601  UT  and  -1620  UT):  this  is  the  effect 
discussed  by  Hall  et  al.  [19911.  But  it  is  also  seen  that  the 
dominant  effect  is  the  minimum  in  Pp  in  the  LLBL:  the 
electrons  make  only  a  minor  contribution  to  the  overall 
pressure  balance.  Also,  the  pressure  contribution  from  high- 
energy  particles  (4(kE<400  keV).  measured  by  the  SULEICA 
instrument  onboard  AMPTE/IRM.  was  negligible  [L.  Kistler, 
private  communication!.  The  occasional  spikes  in  the  P(Qi 
curve  (panel  6)  are  believed  to  be  the  result  of  aliasing 
associated  with  different  sampling  of  field  and  plasma. 

In  summary,  items  to  note  in  this  event  are:  ( 1 )  colocation 
of  field  enhancement  region  and  LLBL;  (2)  gradually 
decreasing  density  of  medium-energy  ions  moving  outward 
from  the  inner  edge  of  the  LLBL;  (3)  dropout  of  medium- 
energy  magnetospheric  electrons  at,  or  earthward  of  the  inner 
edge  of  the  LLBL,  a  dropout  that  is  more  pronounced  than  the 
decrease  in  medium-energy  magnetospheric  ion  densitv  (4) 
near  constancy  of  total  pressure. 

3.  Type  II  Event 

Figure  2  shows  Type  II  behavior,  i.e.,  magnetic  field 
enhancement  in  the  magnetosphere  proper,  for  the  outbound 
pass  of  AMPTE/IRM  on  October  8,  1985.  This  traversal  of 
the  LLBL/MP  region  occurred  near  local  noon  (1120  local 
time)  at  -10.8°  GSE  latitude  and  at  a  geocentric  distance  of 
about  1 1  Re-  At  the  left  edge  of  the  diagram  (0850  UT).  the 
spacecraft  was  in  the  magnetosphere,  proceeding  outward 
toward  the  LLBL.  An  early  gradual  drop  in  electron 
temperature  preceded  the  LLBL  proper  but  other  plasma 
parameters  remained  nearly  constant  there.  However,  a 
gradual  increase  of  field  magnitude  occurred,  from  60  nT  (Pp 
s  1.43  nPa)  at  0853  UT  to  about  70  nT  (Pp  2  1.95  nPa)  when 
the  LLBL  plasma  was  first  encountered,  at  about  0857 : 10  UT. 
At  this  latter  time,  Np  started  to  increase  rapidly,  while  Nje 
dropped  abruptly  and  N2p  more  gradually  to  lower  levels, 
characteristic  of  the  LLBL.  There  was  an  associated  abrupt 
drop  in  Tg  whereas  Tp  decreased  more  gradually  as  Np 
ramped  up  from  the  magnetospheric  to  the  LLBL  level.  At  the 
end  of  this  ramp,  around  0857.45  UT,  the  field  magnitude 
dropped  abruptly;  there  was  little  change  in  the  azimuth  angle, 
cpB,  but  the  elevation  angle,  Xp,  changed  from  essentially  zero 
to  about  +20°. 

One  may  ask  whether  the  spacecraft  entered  the 
magnetosheath  already  at  0857:45  UT,  in  which  case  the 
LLBL  would  have  been  traversed  in  only  about  30  seconds. 
For  a  low-shear  MP  the  possibility  of  misinterpreting 
magnetosheath  magnetic  discontinuities  as  MP  crossings  must 
be  kept  in  mind.  However,  in  the  present  case  there  is 
evidence  to  indicate  that  the  spacecraft  remained  in  the  LLBL 
for  a  long  time,  namely  until  about  0923:40  UT  when  a 
sudden  and  substantial  change  in  field  angle.  tpg,  occurred:  in 
the  interval  0857:45-0923:40  UT,  N2p  and  Nie,  as  well  as  Tp 
and  Te  remained  at  levels  intermediate  between  the 
magnetospheric  and  the  magnetosheath  levels,  a  feature  that  is 
commonly  seen  in  the  LLBL.  Within  10  seconds  after 
0923:40  UT,  N2p  and  Nje  dropped  from  those  intermediate 
levels  to  levels  characteristic  of  the  magnetosheath  with  the 
electron  fluxes  failing  below  the  detection  threshold:  there 
were  associated  decreases  in  Tp  and  Te  as  well  as  changes  in 
the  temperature  anisotropies  or  the  kind  often  seen  at  the  MP. 
If,  as  we  believe,  the  MP  was  traversed  at  0923:40  UT.  it  was 
marked  by  a  substantial  velocity  peak  of  nearly  400  km/s 
compared  to  a  typical  magnetosheath  level  of  about  100  km/s 
and  a  magnetospheric  level  of  50  km/s.  In  fact,  throughout  a 
good  part  of  what  would  then  be  the  LLBL  interval,  from 
0857:45  to  0923:40  UT,  Vp  exceeded  the  flow  speed  in  the 
magnetosheath.  Such  a  long-duration  LLBL  need  not.  and 
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UT  -3.55  GSM  n  =0.977  n  =-.144  n  =-.156  -4.39 

LT  11:02  t  j  i  11:09 

Fig.  2.  AMPTE/IRM  data  during  outbound  pass  through  the 
MP/LLBL  region  on  October  8,  1985;  format  as  in  Fig.  1. 

probably  should  not.  be  interpreted  as  an  indication  of  a  thick 
layer:  there  are  signs  that  the  MP  may  have  reversed  its  radial 
motion  soon  after  0857:45  UT  causing  the  spacecraft  to 
approach  the  inner  edge  of  the  LLBL  around  0901:30  UT. 
There  are  also  indications  of  a  pair  of  MP  encounters  around 
0914  UT  when  a  large  change  in  field  angle,  <pa,  took  place. 
Following  that  change,  a  huge  field  maximum  of  about  85  nT 
is  seen.  It  would  be  difficult  to  incorporate  this  feature  as  a 
semipermanent  pan  of  the  LLBL;  it  seems  more  likely  to  have 
been  a  temporal  effect.  Except  for  this  particular  structure,  it  is 
clear  from  the  figure  that  there  was  a  substantial  diagmagnetic 
field  magnitude  depression,  relative  to  the  magnetospheric 
field,  over  most  of  the  LLBL(0857:45  to  0923:40  UT),  but 
with  a  pronounced  field  maximum  near  its  inner  edge.  This 
field  maximum  and  the  slow  rise  of  the  field,  starting  about 
0853  UT,  and  continuing  until  the  maximum  was  reached  at 
0857:45  UT  is  of  principal  interest  to  us. 

During  the  interval  0853-0857:45  UT,  the  total  pressure 
(panel  6)  gradually  rose  as  a  direct  consequence  of  the  rising 
magnetic  pressure.  In  this  case,  there  is  no  indication  that 
decreases  in  the  thermal  plasma  pressure  were  present  to 
counterbalance  the  increasing  magnetic  pressure:  in  fact,  near 
the  end  of  this  interval,  right  at  the  field  maximum,  the  plasma 
perpendicular  pressure  started  to  increase  instead.  The 
pressure  contribution  from  high-energy  particles  (40<E<400 


keV),  measured  by  SULEICA,  was  not  negligible  for  this 
event:  it  dropped  from  about  0. 1  nPa  prior  to  0850  UT  to 
0.05  nPa  at  0855  UT  and  then  to  zero  at  approximately 
0856:30  UT.  However,  this  decrease  cannot  compensate  for 
the  increase  in  P tot.  shown  in  panel  6  of  Figure  2,  from  about 
1 .6  nPa  at  0853  UT  to  2.4  nPa  at  0857:45  UT.  In  the  regions 
of  depressed  magnetic  field  following  after  the  latter  time,  the 
thermal  plasma  pressure  increased  a  great  deal  (while  the  high- 
energy  particle  pressure  remained  zero)  but  did  not  fully 
compensate  for  the  depressed  magnetic  pressure.  An 
exceptionally  large  maximum  in  P\oi  was  associated  with  the 
aforementioned  large  field  maximum  at  0919:20  UT.  This  is 
an  indication  of  two  or  three-dimensional  structure  and/or 
temporal  evolution  associated  with  this  field  structure. 

The  most  significant  features  of  this  traversal  are:  (1) 
gradual  field  enhancement  in  the  magnetosphere  prior  to  the 
spacecraft  entering  the  LLBL,  with  the  field  maximum  located 
at  the  inner  edge  of  that  layer.  (2)  nearly  constant  level  of  the 
medium-energy  ion  density.  Njp,  in  the  magnetospheric  part 
of  the  field  enhancement  region;  (3)  partial  dropout  of 
medium-energy  magnetospheric  electrons  and  an  associated 
abrupt  drop  in  electron  temperature  near  the  inner  edge  of  the 
LLBL;  (4)  a  smaller  drop  in  N2p  at  the  inner  edge  of  the 
LLBL;  (5)  absence  of  total  pressure  balance  in  the  interval 
0853-0857:45  UT. 

4.  Discussion 

The  near  constancy  of  total  pressure  in  Type  I  MP/LLBL 
crossings  indicates  that  field  curvature  effects  and  temporal 
effects  play  little  role  in  them.  The  principal  item  that  needs  to 
be  examined  is  the  cause  of  the  depressed  densities,  N2p  and 
N2e>  of  medium-energy  magnetospheric  ions  and  electrons  in 
the  LLBL.  It  is  the  depression  in  N2p  that  produces  the  main 
defect  in  total  plasma  pressure  in  the  LLBL,  a  defect  that  is  in 
turn  compensated  for  by  excess  magnetic  pressure  in  that 
layer.  Although  we  agree  with  Neugebauer  et  al.  [  1974]  that 
particles  whose  guiding  centers  are  brought  within  one 
gyroradius  of  the  magnetopause  may  get  lost  to  the 
magnetosheath,  and  although  the  LLBL  could  perhaps  at  times 
be  as  thin  as  a  typical  medium-energy  ion  gyroradius,  we  do 
not  see  how  the  depression  in  the  medium-energy  electron 
density,  N2e<  which  started  substantially  earthward  of  the 
LLBL  in  the  July  3  event  but  which  more  typically  marks  its 
inner  edge,  can  be  accounted  for  in  this  manner.  Possible 
explanations  for  this  behavior  of  N2e  are: 

(a)  Energetic  particle  diffusion  from  the  magnetosphere 
towards  the  MP  could  be  an  important  effect.  From  the 
considerably  steeper  Np  profiles  it  would  appear  that  inward 
diffusion  of  magnetosheath  protons  is  a  much  less  effective 
process  than  outward  diffusion  of  more  energetic  particles. 

(b)  Pitch-angle  scattering  with  associated  particle 
precipitation  could  be  responsible  for  the  depiction  of  N2e  and 
N2p  in  the  LLBL  and  perhaps  for  the  depressed  (relative  to  the 
magnetosheath)  thermal  plasma  density,  Np,  there. 

(c)  The  LLBL  could  be  on  open  field  lines,  i.e.,  field  lines 
with  only  one  end  in  the  ionosphere. 

(d)  Magnetosheath  plasma  could  have  entered  onto  field 
lines  in  the  LLBL  at  some  upstream  location  where  those  field 
lines  were  temporarily  opened  by  reconnection,  allowing 
magnetosheath  plasma  to  enter  onto  them  and  energetic 
magnetospheric  ions  and  electrons  to  be  drained  from  them. 
Subsequently,  the  field  lines  closed  again,  through  a  second 
reconnection  process  (Kan,  1988],  and  were  then  transported 
tailward  along  the  magnetopause  to  the  observation  site.  In  an 
alternate  scenario  [Cowley,  1981],  the  entire  LLBL  is 
incorporated  into  the  magnetosphere  by  cusp  reconnection 
during  periods  of  northward  interplanetary  magnetic  field;  if 
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formed  in  this  manner  the  LLBL  would  be  devoid,  or  partially 
devoid  of  energetic  magnetospheric  particles. 

<e)  Nlagnetosheath  plasma  could  have  moved  onto  closed 
magnetospheric  field  lines  by  E  x  B  drift  in  some  localized 
region  where  the  magnetosheaih  field  was  aligned  with  the 
magnetospheric  field.  In  such  a  process,  the  energetic 
magnetospheric  panicle  population  would  be  pushed  out  of  the 
way.  also  as  a  result  of  E  x  B  drift. 

Without  additional  information,  we  cannot  eliminate  any  of 
these  possibilities.  However,  if  (c),  (d)  or  (e)  were  applicable 
in  their  purest  form.  i.e..  without  the  presence  of  particle 
diffusion,  then  the  density  profiles  N2p  and  N2e  at  the  inner 
edge  of  the  LLBL  would  be  much  steeper  than  they  are 
observed  to  be.  For  this  reason,  we  believe  that  diffusion 
must  play  a  significant  role,  although  it  seems  likely  to  occur 
in  combination  with  one  of  the  other  effects.  Also,  the 
steepness  of  the  Np  profiles  suggests  that  the  diffusion  is  by 
microscopic  rather  than  macroscopic  turbulence,  since  the 
latter  would  operate  equally  as  effectively  on  the  low  energy 
LLBL  plasma  as  on  the  more  energetic  panicles. 

Characteristic  features  of  Type  II  traversals  of  the 
MP/LLBL  region  are  that  the  magnetic-field  magnitude  is 
enhanced  well  within  the  magnetosphere  proper  and  that  the 
total  pressure,  Plot.  is  not  constant  in  this  region.  We  argue 
that  these  effects  are  caused  by  field-line  curvature, 
presumably  associated  with  waves  or  bulges  on  the 
magnetopause  moving  past  the  spacecraft.  It  is  well  known 
that  short-duration  overshoots  in  /’tot  occur  during  flux 
transfer  events  (FTEs).  The  explanation  for  this  effect  is  the 
draping  of  ambient  magnetic-field  lines  around  an  elongated 
MP  structure  such  as  a  flux  tube  (along  with  twisting  of  the 
field  within  the  tube).  Fanrugia  et  al.  (1987]  modeled  the 
draping  effect  in  terms  of  a  vacuum  field,  generated  by  the 
superposition  of  the  uniform  ambient  field  component,  Bqx, 


Fig.  3.  Field-line  aspect  ratio,  lill,  and  maximum  deflection 
angle,  AX,  as  a  function  of  maximum  increase  in  transverse 
field,  ABX.  for  Farrugia  et  al.  [  19H7|  model.  Field  maximum 
occurs  at  A. 

transverse  to  the  flux  tube,  a  constant  field.  Bqj,  along  the 
tube,  and  a  two-dimensional  dipole  field,  representing  the  field 
perturbation  caused  by  the  tube.  This  configuration  has  a  field 
maximum  immediately  above  the  tube.  Using  this  simple 
model,  we  can  estimate  the  field-line  deformation  needed  to 
increase  the  transverse  field  by  an  observed  amount.  The 
result  is  shown  in  Figure  3  in  which  the  ratio  of  height  to 
length,  hll,  of  a  draped  field-line  surface  through  the  field- 


maximum  point.  A,  is  shown  as  a  function  of.  Afit/B<u.  the 
ratio  of  the  enhancement  of  the  transverse  field  component  to 
the  ambient  value  of  that  component.  Also  shown  is  the 
maximum  transverse  field-angle  change,  AX,  observed  bv  a 
spacecraft  as  it  travels  parallel  to  the  magnetopause  towards  the 
field  maximum.  It  is  seen  that  even  gentle  bending  of  the 
field-line  surfaces  may  lead  to  a  substantial  field  increase.  In 
the  October  8,  1985,  event,  a  field  increase  from  60  nT  to  70 
nT  is  observed.  Minimum-variance  analysis  of  the  magnetic- 
field  data  in  the  interval  0853:00-0856:48  UT  indicates  that 
most  of  this  field  was  in  the  transverse  (x)  direction.  Figure  3 
shows  that  such  an  increase  could  be  produced  by  a  curved 
field-line  surface  having  aspect  ratio  li/l  -  1/33.  with  a 
maximum  observed  field  deflection  of  AX  =  6°  occurring  well 
before  the  field  maximum  is  reached.  A  negative  deflection  in 
X/)  of  this  order  of  magnitude  is  in  fact  seen  to  precede  the 
field  maximum  (panel  4  of  Figure  2). 

Even  near  local  noon,  it  is  not  surprising  to  find  the  MP  to 
experience  frequent  undulations  of  this  magnitude:  multiple 
encounters  of  a  spacecraft  with  the  LLBL/MP  would  be  a 
natural  consequence.  Thus  Type  II  behavior  should,  and 
indeed  frequently  is  observed  to,  accompany  such  multiple 
encounters  (e.g..  Figure  9  of  Neugebauer  et  al.  (1974)). 
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